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ABSTRACT 

The  point  charge  model  is  used  to  calculate  the  crystal  structure  of  sulfuric 
acid  (H2SO4)  with  the  6-3 IG**  basis  set.  The  point  charge  model  accurately 
reproduces  the  structural  trends  which  occur  in  transforming  from  the  gas  to  the 
solid  phase.  The  calculated  crystal  structure  of  sulfuric  acid  is  in  reasonably  good 
agreement  with  both  the  X-ray  and  neutron-diffraction  structures. 

The  point  charge  model  is  shown  to  precisely  simulate  the  deformation 
forces  which  are  present  in  the  solid  upon  crystallization.  The  point  charge  model 
exhibits  a  definite  shift  of  electron  density  from  the  bridging  hydrogens  to  the 
acceptor  atoms,  identical  to  those  found  in  other  ab  initio  studies.  The  calculated 
crystal  structures  are  insensitive  to  the  magnitude  of  the  point  charges. 

Two  new  iterative  techniques,  using  the  point  charge  model,  are  introduced 
which  give  superior  results  to  any  of  the  single  optimization  cycle  methods.  These 
iterative  techniques  account  for  any  forces  which  are  not  electrostatic  by  nature. 

The  gas  phase  structure  of  H2SO4  is  optimized  using  the  ST0-3G*,  4- 
31G**,  and  6-3 IG**  basis  sets.  Comparisons  are  made  between  the  programs 
TEXAS,  TEXAS(D),  and  Gaussian-82. 
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CHAPTER  ONE 

INTRODUCTION 

Over  the  years,  many  researchers  have  shown  the  ability  to  apply  quantum 
theory  to  chemistry  in  order  to  supplement  experimental  results.  However,  it  has 
only  been  in  the  past  decade  with  the  great  advances  in  computers  coupled  with  the 
development  of  efficient  algorithms  that  ab  initio  calculations  have  become  practical 
on  molecules  of  more  than  a  few  atoms. 

The  term  "ab  initio"  is  frequently  cited  in  the  literature  but  rarely  defined. 
Ab  initio  calculations  can  basically  be  thought  of  as  those  which  treat  all  the 
electrons  and  do  not  utilize  parameters  which  are  fit  to  experimental  data  as 
compared  to  semi-empirical  calculations.  For  example,  a  Hartree-Fock  calculation 
of  a  restricted  closed-shell  single  determinant  wave  function  is  one  in  which  no 
approximation  is  made  to  the  integrals  or  the  electronic  Hamiltonian.  However,  it 
is  completely  specified  by  the  choice  of  a  basis  set  and  the  coordinates  of  the 
nuclei.  1 

One  area  in  which  ab  initio  calculations  have  been  particularly  successful  is 
in  evaluation  of  molecular  structure.  The  theoretical  determination  of  a  molecular 
structure  involves  minimizing  the  total  energy  of  the  molecule  with  respect  to 
simultaneous  variation  of  the  intemal  coordinates.^  The  application  of  this 
procedure  to  complex  molecules  became  practical  only  after  the  introduction  of  the 
gradient  method  by  Pulay.  This  method  gready  reduced  the  time  to  calculate  an 
optimum  geometry  especially  for  large  molecules.^  The  resulting  optimized 
structure  has  a  known  basis  set  dependance  which  results  in  a  constant  error  for  a 
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given  parameter  called  the  offset  value.  This  offset  value  is  reasonably  constant 
over  a  wide  range  of  molecules.'^  Boggs  and  coworkers  have  reported  extremely 
accurate  structures  by  combining  experimental  and  theoretical  data.  The  structures 
presented  in  their  study  are  more  accurate  than  either  theoretical  or  experimental 
techniques  could  produce  alone.^  Not  only  have  ab  initio  calculations  been  able  to 
accurately  reproduce  experimental  bond  lengths,  bond  angles,  and  conformations 
but  in  some  cases  ab  initio  calculations  have  actually  shown  where  experiments 
have  determined  incorrect  structural  parameters.  Although  the  latter  is  the 
exception,  it  shows  both  the  value  and  validity  of  ab  initio  calculations  as  a 
supplement  to  experimental  data. 

The  ability  to  determine  accurate  molecular  structures  is  the  key  to  many 
areas  of  chemistry.  The  structure  of  a  molecule  not  only  defines  its  physical 
properties  but  also  influences  how  it  reacts  with  other  molecules.  This  information 
is  vital  to  such  diverse  areas  as  the  development  of  "stable",  powerful  explosives 
and  propellants,  the  synthesis  of  effective  catalysts,  the  understanding  of  how 
enzymes  with  only  specific  structures  can  fit  into  certain  substrates,  the 
determination  of  complex  reaction  mechanisms,  and  the  production  of 
superconductors.  However,  to  a  chemist,  the  most  fundamental  concern  which 
structural  information  can  shed  Ught  on  is  chemical  bonding.^  The  understanding 
of  chemical  bonding  is  fundamental  to  all  areas  of  science. 

One  area  in  which  there  is  a  great  deal  of  interest  is  the  structural  differences 
which  occur  in  some  molecules.  These  differences  can  be  between  the  same  type 
of  bond  (eg,  C-H)  or  functional  group  (eg,  CH3)  over  a  variety  of  similar 
compounds  with  different  substituents.  These  changes  can  also  occur  within  the 


same  molecule  when  it  is  in  a  different  phase  or  chemical  environment  (eg,  in 
solution).  These  structural  differences  need  not  be  large  to  be  important.  In  fact, 
they  can  be  very  small.  Differences  on  the  order  of  0.01  A  in  a  bond  length  or  a 
couple  of  degrees  in  a  bond  angle  may  reveal  extremely  valuable  structural 
information  if  they  are  reliably  determined.  However,  only  real  structural 
differences  should  be  considered  and  any  apparent  differences  which  are  detected 
must  be  eliminated  if  meaningful  comparisons  are  to  be  made.^  Apparent 
differences  can  arise  from  either  the  particular  interaction  of  radiation  with  matter  or 
because  of  the  vibrational  averaging  which  occurs  in  a  particular  experiment.  It  is 
important  to  note  that  all  experimental  techniques  produce  geometrical  parameters 
which  are  based  on  averaging  of  the  molecular  vibrations.  However,  at  least 
theoretically,  ab  initio  calculations  give  equilibrium  intermolecular  distances  which 
correspond  to  the  minimum  of  the  potential  energy  surface  and  are  not  subject  to 
molecular  vibrations. 

The  vast  majority  of  all  ab  initio  calculations  are  performed  on  isolated 
"gas"  phase  molecules.  The  reason  for  this  is  quite  simple  -  computer  capacity. 
The  most  accurate  results  come  from  ab  initio  calculations  performed  with  large 
basis  sets.  However,  current  computer  capacity  limits  large  basis  set  calculations  to 
no  more  than  about  two  dozen  atoms.  This  limit  of  two  dozen  atoms  with  a  fairly 
large  basis  set  is  pushing  the  capacity  of  even  the  Cray  supercomputer.  Thus,  in 
order  to  conduct  conventional  ab  initio  calculations  on  solids,  the  size  of  the  basis 
set  would  have  to  be  greatly  reduced,  a  lower  level  of  theory  must  be  used,  or 
assumptions  must  be  made  about  the  symmetry  of  the  molecule.  These  restrictions 
almost  always  result  in  less  accurate  molecular  parameters. 


The  intent  of  my  work  was  to  study  the  structural  changes  which  occur  in  a 
molecule  between  the  gas  and  soUd  phase  through  the  use  of  ab  initio  calculations. 
In  selecting  a  particular  molecule  to  study,  I  estabUshed  several  criteria.  First,  the 
number  of  atoms  in  the  molecule  must  be  small  enough  that  a  large  basis  set  could 
be  used  to  give  accurate  results  but  be  within  the  capacity  of  the  Cray 
supercomputer.  All  of  the  calculations  in  this  study  were  performed  on  the  Cray  X- 
MP/24  at  the  University  of  Texas  Center  for  High  Performance  Computing 
(UTCHPC).  Second,  there  must  be  an  accurate  experimental  determination  of  the 
gas  and  solid  phase  structure  so  that  comparisons  could  be  made  with  the 
calculations.  Finally,  if  possible  the  molecule  should  be  one  which  is  of  general 
interest.  After  careful  consideration,  sulfuric  acid  (H2SO4)  was  selected  as  the 
molecule  to  study  which  met  all  of  the  above  criteria. 

In  order  to  reduce  the  calculation  of  solids  to  one  that  can  be  handled  by 
present  computers,  a  simpler  method  must  be  used  in  order  to  simulate  the  effects 
of  the  crystal  lattice.  The  problem  is  basically  how  to  accurately  represent  the 
molecular-charge  distribution  of  the  surrounding  molecules  in  the  solid.  The 
simplest  is  to  use  the  Mulliken  population  analysis.^  The  Mulliken  population 
analysis  comes  from  parameters  that  are  already  calculated  in  quantum  mechanical 
calculations.  The  major  problem  with  this  method  is  that  the  assignment  of  charges 
to  the  atoms  is  rather  arbitrary  and  is  often  very  basis  set  dependent. ^  However,  a 
previous  study  done  using  the  Mulliken  population  analysis  and  point  charges  was 
able  to  accurately  reproduce  the  crystal  structure  of  cyanoformamide 
(NCC0NH2).^  The  purpose  of  my  study  was  to  apply  the  point  charge  method  to 
a  different  type  of  molecule  and  to  develop  some  new  techniques  in  order  to  see 


their  effect.  The  theoretical  methods  used  in  this  study  are  presented  in  Chapter  2. 
Chapter  2  includes  not  only  the  basis  set  and  level  of  theory  considerations  but  also 
a  review  of  the  various  methods  which  have  been  used  to  simulate  crystal  field 
effects.  Chapter  3  contains  the  sulfuric  acid  results  including  the  optimized  gas 
phase  geometry  and  the  various  techniques  used  to  calculate  the  solid  geometry. 
Chapter  4  presents  my  conclusions  and  suggestions  for  future  studies. 


CHAPTER  TWO 

METHODS 

In  any  ab  initio  calculation,  the  selection  of  the  basis  set  and  the  level  of 
theory  are  crucial.  A  decision  must  be  made  between  the  computational  cost  and 
the  desired  accuracy.  Obviously,  the  most  accurate  results  are  desired  from  any 
calculation.  However,  a  structure  accurate  to  0.01  A  at  one-hundredth  the  cost  of 
one  which  is  accurate  to  0.001  A  maybe  preferred  if  the  fu^st  structure  demonstrates 
the  effects  which  are  being  studied.  Although,  this  is  a  hypothetical  example,  it 
points  out  the  importance  of  extremely  careful  selection  of  the  basis  set  and  level  of 
theory.  In  my  particular  study,  the  selection  of  a  method  to  simulate  the  crystal- 
field  effects  is  also  a  major  factor.  Again,  selection  involves  a  compromise.  In  this 
case,  a  compromise  must  be  made  between  computational  cost,  accuracy,  and  ease 
of  implementation.  In  this  chapter,  I  will  discuss  the  various  methods  used  in  this 
suidy.  The  level  of  theory  and  basis  set  are  discussed  first,  followed  by  the  crystal- 
field  simulation  method.  Finally,  I  will  discuss  the  convergence  criteria  used  in  this 
study. 

LEVEL  OF  THEORY  AND  BASIS  SET  SELECTION 

The  ultimate  goal  of  quantum  chemistry  is  the  calculation  of  the  solution  to 
the  Schrodinger  equation  resulting  in  an  exact  molecular  wave  function.  However, 
that  requires  an  infinitely  large  basis  set  expansion  of  the  wave  function  along  with 
the  implementation  of  full  configuration  interaction.  This  is  clearly  an  unattainable 
goal.  In  actuality,  the  basis  set  is  truncated  to  a  finite  size  and  the  electron 
correlation,  if  treated  at  all,  is  hmited  to  a  few  configurations.  These  specifications 
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define  a  theoretical  model  within  which  all  structures,  energies,  and  other  physical 
properties  can  be  studied.  ^^  While  proper  selection  of  the  model  is  important, 
recognizing  and  understanding  the  limitations  of  the  given  model  on  the  results  is 
imperative. 

The  main  sources  of  error  in  any  geometry  optimization  using  the  gradient 
method  are  neglect  of  electron  correlation  (  where  it  is  either  ignored  or 
incompletely  treated)  and  the  use  of  a  finite  basis  set.  The  combined  effect  of  these 
errors  on  the  results  obtained  from  a  computed  molecular  structure  are  shown  in 
Figure  2- 1 .  The  figure  depicts  the  error  in  some  arbitrary  structural  parameter 
(such  as  a  C-H  bond)  in  a  variety  of  molecular  environments  plotted  against 
increasing  basis  set  size.  When  an  infinite  basis  set  is  used  (the  Hartree-Fock 
limit),  the  error  due  to  the  neglect  of  electron  correlation  is  found  to  be  remarkably 
constant  for  the  given  parameter  over  a  wide  range  of  molecules.  For  much  smaller 
basis  sets,  the  error  is  expected  to  fall  within  the  shaded  area.  However,  at  some 
point  (X),  the  basis  set  is  sufficiently  large  that  the  scatter  in  the  calculated 
parameter  is  extremely  small.  This  results  in  a  constant  known  error  for  the  given 
parameter  called  the  offset  value.  The  remaining  absolute  error  is  due  to  the 
truncation  of  the  basis  set  and  the  neglect  of  electron  correlation.  Thus,  the  major 
flaw  with  small  basis  sets  is  the  inconsistency  of  the  error.'* 

In  order  to  keep  the  computational  costs  low  for  large  molecules,  the 
calculation  should  be  performed  as  close  as  possible  to  the  point  X.  Any  increase 
in  the  basis  set  results  in  an  increase  in  computation  time  and  disk  storage  space 
proportional  to  the  fourth  power  of  the  number  of  basis  functions  (N^).  Beyond 
point  X,  the  wave  function  approaches  the  Hartree-Fock  limit  where  only  electron 
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correlation  error  remains.  The  incorporation  of  any  high-level  electron  correlation 
treatments  results  in  a  large  increase  in  computation  time  and  disk  storage  roughly 
proportional  to  n^  where  n  is  the  number  of  electrons  in  the  molecule.  ^  ^  Even  the 
most  efficient  correlation  methods,  such  as  the  local  correlation  method  ^  2,  are  two 
to  three  times  more  expensive  than  a  single  determinant  calculation.  Therefore, 
depending  upon  the  size  of  the  molecule  and  the  number  of  basis  functions,  the 
treatment  of  electron  correlation  may  be  prohibitively  expensive  for  a  given 
molecule. 

Based  on  the  above  observation,  calculations  can  be  done  efficiently  at  point 
X.  For  many  structural  parameters,  point  X  corresponds  to  approximately  double 
zeta  basis  sets.  However,  for  some  parameters  such  as  torsional  angles  around 
nitrogen  or  oxygen,  it  is  essential  to  add  at  least  one  set  of  polarization  functions  to 
the  basis  to  reach  that  point.  The  offset  value  for  most  bond  angles  is  zero  using 
this  basis  set.  However,  bond  distances  have  non-zero  offset  values  at  this  level 
and  for  some  types  of  bonds  the  basis  set  must  go  beyond  double  zeta.^^  \ 
common  misconception  is  that  the  addition  of  d  orbitals  (polarization  functions)  to 
first  and  second  row  elements  is  incorrect  since  d  orbitals  in  these  elements  are  not 
actually  occupied.  The  addition  of  polarization  functions  is  only  done  to  allow 
greater  flexibility  in  the  description  of  the  bonds.  Their  addition  does  not  imply  that 
d  orbitals  are  actually  involved  in  the  formation  of  these  bonds. 

Because  of  the  size  of  the  H2SO4  molecule,  the  treatment  of  electron 
correlation  is  prohibitively  expensive.  Unless  stated  otherwise,  the  calculations  in 
my  study  were  completed  at  the  single  determinant  level  using  the  gradient  program 
TEXAS  1'*  with  no  symmetry  restrictions.  During  the  course  of  my  study,  Pulay's 
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group  completed  a  new  version  of  TEXAS  which  was  used  for  a  few  additional 
calculations.  The  major  improvements  are  increased  speed  and  the  ability  to  use 
true  Gaussian  d  functions.  After  some  optimization  of  the  program  by  Cordell  for 
the  Cray,  I  have  seen  up  to  a  30%  increase  in  speed  over  the  old  version  of 
TEXAS.  The  inclusion  of  true  d  functions  has  replaced  the  use  of  displaced  p 
functions  which  were  used  to  construct  d  functions  in  the  previous  version  of 
TEXAS.  Throughout  this  paper,  the  new  version  of  TEXAS  will  be  referred  to  as 
TEXAS  (D). 

As  stated  earlier,  one  of  the  most  important  parts  of  any  ab  initio  calculation 
is  the  selection  of  a  basis  set.  The  literature  is  inundated  with  almost  100  basis 
sets^5  which  range  from  the  minimal  ST0-3G  basis  set  to  very  extensive  triple  zeta 
plus  multiple  polarization  function  basis  sets.  Thus,  the  selection  of  a  basis  set  for 
a  particular  problem  must  be  made  based  on  the  ab  initio  program  used,  the 
accuracy  required,  and  how  much  computation  time  is  available.  The  TEXAS 
program  efficientiy  uses  split  valence  shell  basis  sets  where  the  s  and  p  Gaussian 
exponents  are  equal.  Pople's  6-3 IG**  basis  sets  are  used  almost  exclusively 
throughout  this  study.  16,17, 18,19,20  Some  calculations,  which  will  be  specifically 
mentioned,  were  completed  using  a  4-3 IG**  basis  set  on  sulfur^^  with  a  6-3 IG** 
basis  set  on  oxygen  and  hydrogen  (referred  to  throughout  the  remainder  of  this 
paper  as  4-3 IG**).  The  asterisks  in  the  notation  for  the  basis  sets  refer  to  the 
inclusion  of  polarization  functions.  The  first  asterisk  refers  to  the  addition  of  d 
functions  to  the  heavy  atoms  sulfur  (exponent=0.65)'*7  and  oxygen 
(exponent=0.8)22.  The  second  asterisk  refers  to  the  addition  of  p  functions  to 
hydrogen  (exponent=  1.1)22.  As  discussed  earlier,  the  inclusion  of  polarization 
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functions  is  essential  to  accurately  obtain  torsional  angles  around  oxygen.  Another 
basis  set  which  was  used  to  study  the  sulfuric  acid  dimer  was  Pople's  STO-SG*.^^ 
That  basis  set  consists  of  the  minimal  STO-3G  basis  with  the  addition  of  five  d 
functions  only  to  second  row  atoms. 

METHOD  TO  SIMULATE  CRYSTAL-FIELD 

Numerous  studies  have  been  made  in  order  to  apply  ab  initio  calculations  to 
solids.  One  reason  is  the  advantages  which  ab  initio  calculations  offer  over  some 
experimental  methods.  While  X-ray  and  neutron  diffraction  provide  valuable 
information  about  the  spatial  arrangement  of  the  atoms  in  molecules  and  the  packing 
of  the  molecules  in  lattices,  they  tell  us  nothing  of  the  forces  which  determine  the 
crystal  structure.  In  addition,  they  give  no  information  on  why  the  molecules  pack 
in  the  observed  space  group  or  the  influence  of  crystal  forces  on  the  conformation 
of  flexible  molecules.^^  In  some  cases  the  experimental  methods  cannot  "locate" 
particular  parts  of  a  molecule,  especially  hydrogens.  In  the  case  of  sulfuric  acid, 
three  X-ray  diffraction  studies,  with  the  latest  being  done  in  1978,  could  not 
"locate"  the  hydrogens.  It  was  not  until  a  neutron  diffraction  study  was  done  in 
1983  that  the  hydrogens  were  located.  Another  problem  with  the  experimental  data 
is  the  vibrational  averaging  which  is  inherent  to  all  experimental  methods.  It  is 
these  questions  and  others  to  which  ab  initio  calculations  can  be  directly  applied. 

One  area  in  which  theoretical  methods  have  been  frequentiy  applied  is  the 
study  of  hydrogen-bonding.  The  ab  initio  methods  which  have  been  used  range 
from  the  complex  use  of  an  atom-centered  multipole  expansion's  to  the  simple 
point  charge  model.^  It  was  not  until  the  study  by  Saeb0,  et  al.  that  the  influence  of 
crystal  forces  on  the  intramolecular  geometry  was  studied  in  detail.  Saeb0  reported 
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impressive  agreement  between  X-ray  diffraction  and  calculated  stmcuiral 
parameters  by  using  point  charges  based  on  a  Mulliken  population  analysis. 
Regardless  of  the  methods  which  have  been  used,  the  importance  of  the 
electrostatic  interaction  energy  in  determining  the  energetically  most  favorable 
geometry  is  well  established.^'^'^'^ -26,27  j^  fact,  it  led  Bonaccorsi,  et  al  28  to  the 
"electrostatic  assumption".  Bonaccorsi  and  coworkers  postulated  that  the  minima 
in  the  total  interaction  energy  coincide  with  those  in  the  electrostatic  energy  for 
rotational  degrees  of  freedom.28  Thus,  it  is  clear  that  approximating  the 
electrostatic  energy  via  a  model  can  reveal  valuable  molecular  structure  information 
if  the  model  is  sufficientiy  detailed. 

One  study  which  specifically  considered  crystal-field  effects  was  the 
calculation  of  the  rotational  conformation  of  acetamide.29  In  that  case,  the  crystal 
structure  had  shown  that  the  torsion  angle  was  between  0°  and  30°.  However, 
several  attempts  by  other  authors  to  compute  the  lattice  energies  and  total 
crystallographic  conformational  energies  (lattice+torsional)  could  not  account  for 
that  torsional  angle.  A  procedure  was  developed  to  compute  the  lattice  energy  as 
the  sum  of  three  long-range  contributions  (electrostatic,  polarization,  and 
dispersion)  and  a  short-range  repulsive  contribution.  That  procedure  yielded  a 
minimum  lattice  energy  with  a  torsional  angle  between  0°  and  30°  as  observed  in  the 
crystal  structure.29  Similarly,  through  the  use  of  calculated  electrostatic  molecular 
potentials,  researchers  were  able  to  compute  the  preferred  arrangement  of  linear 
chains  of  formamide  molecules  as  observed  in  the  crystal.^O 

The  combination  of  experimental  and  theoretical  methods  has  been  used  in 
several  studies.  One  study  of  hydrogen  bonding,  crystal  packing,  and  the  effect  of 
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crystal  forces  on  the  molecular  conformation  of  amides  and  carboxylic  acids 
expounded  upon  the  "marriage"  of  X-ray  diffraction  and  ab  initio  calculations  24 
The  goal  of  that  study  was  to  obtain  a  general  force  field  for  organic  and  biological 
molecules.  In  that  study,  the  lattice  energy  was  calculated  as  a  sum  of  the  Lennard- 
Jones  potential  and  the  Coulombic  interaction.  The  researchers  were  able  to  use 
their  procedure  to  explain  the  crystal-packing  mode  of  formic  and  acetic  acid.  Their 
study  clearly  showed  that  the  electrostatic  interactions  were  the  overriding  factor  in 
determining  the  particular  packing  occurring.24  In  another  "marriage", 
experimental  electron  densities  were  combined  with  calculations  in  order  to  study 
the  role  of  Coulomb  forces  in  the  crystal  packing  of  amides.  In  this  case,  the 
researchers  were  able  to  use  the  Coulomb  interaction  energies  in  order  to  determine 
the  possible  and  preferred  molecular  packing  modes.  In  addition,  they  were  able  to 
obtain  reliable  estimates  of  the  van  der  Waals  atom-atom  potential  parameters  using 
the  calculated  Coulomb  energies.^^ 

In  one  investigation,  the  electrostatic  energy  in  the  formic  acid  crystal  was 
studied.26  The  geometry  of  formic  acid  was  not  optimized;  instead  the  researchers 
used  the  experimental  gas  and  crystal  geometries  to  conduct  their  study.  In  their 
study,  two  approaches  were  pursued  -  multipole  and  point  charge  lattice  sums.  In 
the  multipole  method,  the  electrostatic  energy  is  expanded  in  a  series  to  yield  the 
various  multipole  moments.  That  particular  study  used  multipole  moments  up  to 
the  sixth  moment.  The  electrostatic  energy  then  becomes  a  sum  of  the  interaction 
energy  between  the  polar  multipole  moments.  In  the  point  charge  method,  the 
charge  distribution  is  represented  by  a  set  of  point  charges  which  were  obtained 
from  a  MuUiken  population  analysis.  The  electrostatic  energy  simply  reduces  to  the 
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sum  of  the  interaction  energy  between  the  point  charges.  Several  important 
conclusions  were  made  from  the  results.  First,  the  electrostatic  portion  of  the  lattice 
energy  of  formic  acid  was  both  large  in  magnitude  and  varied  strongly  with  the 
orientation  of  the  molecule.  Several  orientations,  including  the  one  observed 
experimentally,  were  found  to  minimize  the  electrostatic  energy.  The  location  of 
these  minima  were  dictated  by  interactions  with  nearest-neighbor  molecules. 
Second,  the  point  charge  method  agreed  qualitatively  with  the  multipole  series  for 
the  dependence  of  the  electrostatic  energy  on  the  molecular  orientation.  Finally  and 
most  importantiy,  the  electrostatic  component  of  the  lattice  energy  was  found  to  act 
as  the  driving  force  for  molecular  distortion  upon  crystallization.26 

The  conclusion  which  can  be  drawn  from  all  of  the  above  studies  is  that 
theoretical  calculations  can  give  extremely  valuable  information  about  a  wide  variety 
of  crystal  parameters.  However,  it  is  now  time  to  specifically  address  the  crystal 
parameter  which  is  of  interest  in  my  study  -  molecular  structure.  A  straightforward 
approach  which  has  been  appUed  to  the  calculation  of  molecular  structure  is  to 
compute  it  directiy  in  the  sohd.  This  involves  selecting  a  model  of  an  extended 
"supermolecule"  that  is  sufficientiy  large  to  represent  the  condition  existing  inside 
the  crystal.  That  same  method  has  been  previously  used  to  study  such  phenomena 
as  the  adsorption  of  a  gas  on  a  surface,  where  the  solid  is  represented  by  a  small 
number  of  atoms.  One  way  of  implementing  this  method  is  to  consider  all  the 
molecules  within  a  unit  cell  as  the  "supermolecule"  and  then  impose  boundary 
conditions  to  correct  for  effects  at  the  cell  edge.  However,  as  the  size  of  the 
molecular  system  increases,  the  accuracy  of  the  computed  results  suffers.  As  a 
result,  the  "supermolecule"  method  has  not  been  conducive  to  highly  accurate 
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results  for  molecules  of  appreciable  size.^  A  recent  study  using  that  method  has 
been  completed  by  Cordell  and  Boggs  on  furan  molecules.^  1  The  calculations  were 
done  on  furan  dimers  and  trimers  located  in  the  spatial  orientation  found  in  the 
crystal.  Their  study  was  able  to  show  which  neighbor  molecules  stabilized  the 
crystal  structure  and  to  what  extent.  The  interaction  energies  of  the  various  subsets 
of  the  unit  cell  were  shown  to  be  essentially  additive,  with  the  packing  forces 
expected  to  decrease  the  calculated  crystal  stabiUzation  by  5  to  10%.  In  addition, 
they  concluded  that  the  short  C=C  bond  length  reported  for  the  crystal  was  an 
artifact  of  the  X-ray  diffraction  experiment.^!  However,  as  stated  earlier,  the 
"supermolecule"  technique  is  only  useful  for  relatively  small  molecules.  There  is 
also  an  appreciable  cost  associated  with  it  which  increases  as  N^.  Thus,  large  basis 
set  and/or  large  molecule  calculations  using  the  "supermolecule"  concept  are  very 
expensive. 

Based  on  the  previous  discussion,  there  are  major  problems  which  afflict 
many  of  the  methods  which  have  been  used  to  simulate  crystal-field  effects.  All  of 
the  above  methods  suffer  from  one  or  more  of  the  following  problems:  high 
computational  cost,  low  accuracy,  or  difficulty  in  implementation.  Recently,  a 
method  has  been  proposed  by  Saeb0,  et  al.  to  calculate  crystal  structures  which 
seems  to  overcome  these  problems.  Their  study  consisted  of  investigating  the 
geometry  changes  which  occur  in  cyanoformamide  in  going  from  the  gaseous  to  the 
solid  states,  with  emphasis  on  the  effect  of  the  crystal  forces  on  the  intramolecular 
geometry.  Their  approach  consisted  of  using  a  simple  point  charge  model  to 
simulate  the  crystal-field  effects.  First,  the  geometry  of  the  free  molecule  was 
optimized  and  a  wave  function  obtained.  Next,  a  MuUiken  population  analysis  was 
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nm  to  generate  the  magnitude  of  the  point  charges.  The  point  charges  were  then 
placed  around  the  free  molecule  in  the  spatial  arrangement  which  corresponded  to 
the  previously  determined  X-ray  diffraction  structure.  Since  the  hydrogens  were 
not  located  in  the  X-ray  diffraction  experiment,  the  N-H  bond  distances  were  set  to 
1.0  A.  Finally,  the  geometry  of  the  free  molecule  was  optimized  inside  of  the  point 
charges.  The  resulting  optimized  structure  yielded  bond  lengths  and  bond  angles 
with  a  maximum  deviation  of  0.01 1  A  and  0.4°,  respectively,  compared  to  the  X- 
ray  diffraction  structure.  The  major  advantages  of  their  method  were  not  only  the 
accurate  results  and  the  ease  of  implementation  but  also  the  very  slight  increase  in 
computation  time  compared  to  the  free  molecule.  The  small  increase  in  computation 
time  was  only  due  to  the  increase  in  the  number  of  one-electron  integrals  and 
nuclear  repulsion  terms.9  Thus,  their  method  appears  to  possess  all  the  advantages 
of  the  previously  discussed  methods  with  none  of  the  disadvantages. 

There  are  a  few  problems  with  the  point  charge  model  which  need  to  be 
considered.  The  model  assumes  that  the  crystal-field  effects  can  be  simulated 
through  purely  electrostatic  interactions.  However,  these  are  not  the  only  forces 
which  bind  molecules  together  in  solids.  In  order  to  fully  understand  the  types  of 
forces  which  bind  crystals  together,  the  various  crystals  which  exist  must  be 
reviewed.  There  are  basically  four  types  of  crystals  -  ionic,  covalent,  metallic,  and 
molecular.  Metals  represent  a  special  class  of  solids  and  will  not  be  discussed 
further.  In  ionic  crystals,  the  atoms  are  bound  together  by  largely  electrostatic 
forces  and  there  is  no  such  thing  as  an  individual  molecule.  Covalent  solids  are 
actually  one  huge  molecule  which  is  composed  of  covalently  bonded  atoms.    The 
atoms  are  held  together  by  sharing  electrons  with  their  neighbors.  On  the  other 
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hand,  molecular  crystals  actually  contain  a  collection  of  distinguishable  molecules 
which  are  bound  together.  The  molecules  are  held  together  mostly  by  dipole 
interactions,  van  der  Waals'  forces  and  hydrogen  bonds.  In  reality,  very  few 
crystals  fall  stricdy  into  one  of  these  categories.  Instead,  many  crystals  are  held 
together  by  more  than  one  of  these  forces.  Not  only  are  there  numerous  examples 
of  crystals  composed  of  both  ionic  and  covalent  bonds  but  even  molecular  crystals 
may  contain  a  small  amount  of  ionic  or  covalent  bonding.32  But,  all  of  these  forces 
share  a  common  denominator  -  they  all  depend  upon  charged  particle  interactions. 
In  summary,  the  forces  which  bind  crystals  together  are  very  complex  and  in 
certain  cases  may  not  be  fully  described  by  only  electrostatic  interactions. 

The  foremost  problem  with  the  point  charge  model  used  by  Saeb0  and 
coworkers  is  that  it  does  not  specifically  allow  for  intermolecular  electron  transfer. 
The  neglect  of  electron  transfer  in  their  method  coiJd  certainly  degrade  the  accuracy 
of  calculated  crystal  structures  in  which  covalent  bonding  is  a  major  force. 
However,  since  few  covalent  crystals  exist,  their  method  should  give  satisfactory 
results  for  a  wide  variety  of  other  crystals.  In  fact,  the  study  by  Saeb0  indicated 
that  hydrogen  bonding  was  dominated  by  electrostatic  effects.  WhUe  the  neglect  of 
electron  transfer  must  not  be  forgotten  when  interpreting  results  using  this  method, 
the  point  charge  method  should  produce  reasonable  results  when  applied  to  a  large 
variety  of  crystals,  especially  molecular  crystals  involving  hydrogen  bonding. 

There  still  remain  a  couple  of  problems  with  the  Mulliken  population 
analysis  which  need  to  be  considered.  First,  it  tends  to  produce  charges  which  are 
rather  arbitrary  and  basis  set  dependent.  Any  defects  in  the  point  charge  model  are 
dominated  by  possible  errors  in  the  magnitude  of  the  point  charges.^  Second,  since 
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it  is  based  on  a  simplified  model  of  describing  the  electron  distribution,  the 
Mulliken  population  analysis  often  yields  multipole  moments  which  are  quite 
different  than  those  calculated  from  the  actual  wave  function.^  Although,  there  are 
other  standard  population  analysis  methods  such  as  that  of  Lowdin,  they  suffer 
from  the  same  disadvantages.  ^ 

In  order  to  obtain  more  accurate  charges,  a  different  method  has  been 
recentiy  proposed.^  In  that  method,  the  molecular  electrostatic  potential  evaluated 
at  points  in  space  around  the  molecule  is  used  as  a  guide  and  is  then  fit  to  the  point 
charge  models.  The  procedure  provides  for  selected  restraints  on  the  calculated 
dipole  or  quadrupole  moments,  and  for  the  use  of  additional  point  charges  to 
represent  lone  pairs  in  the  molecule.  That  study  found  very  good  agreement  with 
experimental  enthalpies  of  formation  for  hydrogen  bonding  on  nucleic  acid  pairs. 
Based  on  the  above  results,  the  researchers  found  their  method  superior  to  the 
Mulliken  population  analysis.  Even  given  a  set  of  point  charges  which  accurately 
reproduce  the  molecular  electrostatic  potentials,  it  is  not  a  trivial  task  to  determine  a 
set  of  nonbonded  parameters  that  can  be  used  with  the  charges  to  give  good 
geometries.8  This  method  also  tends  to  assign  larger  absolute  values  of  charge  to 
the  atoms  compared  to  the  Mulliken  population  analysis.  However,  ease  of 
implementation  was  the  overriding  factor  in  choosing  between  these  two  methods. 
While  their  method  may  produce  a  more  accurate  point  charge  model  than  the 
Mulliken  method,  it  is  by  no  means  a  trivial  task  to  implement  their  method.  The 
method  which  provides  the  best  compromise  based  on  the  selection  criteria  is  the 
point  charge  model  using  the  Mulliken  population  analysis. 

Both  versions  of  TEXAS  have  been  modified  to  accommodate  the  point 
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charge  model.  The  programs  require  that  the  point  charges  be  input  in  Cartesian 
coordinates  with  their  associated  net  charge.  The  energy  which  is  calculated 
excludes  the  point  charge-point  charge  interactions  but  includes  the  point  charge- 
molecule  interactions.  Therefore,  the  calculated  energy  accurately  shows  the 
stabilization  of  the  crystal  when  compared  to  the  free  molecule  energy. 

The  major  reason  that  the  Mulliken  population  analysis  is  easy  to  implement 
is  because  the  required  parameters  are  already  calculated  in  any  standard  SCF 
program.  The  only  required  parameters  are  the  atomic  number,  the  density  matrix 
(P),  and  the  overlap  matrix  (S).  The  number  of  electrons  that  are  associated  with 
any  basis  function  is  simply  the  diagonal  element  of  the  product  of  the  density  and 
overlap  matrices  which  corresponds  to  that  basis  function.  The  corresponding 
number  of  electrons  that  are  associated  with  a  given  atom  in  a  molecule  is  obtained 
by  summing  over  all  basis  functions  on  that  atom.  This  assumes  that  all  the  basis 
functions  are  centered  on  atomic  nuclei.  The  net  charge  of  any  atom  (point  charge 
magnitude)  is  simply  the  difference  between  the  charge  of  the  atomic  nucleus  and 
the  number  of  electrons  associated  with  that  atom.  The  net  charge  is  positive  if  the 
number  of  electrons  associated  with  the  atom  is  less  than  the  nuclear  charge  and 
negative  otherwise.  Although  this  method  is  simple,  it  must  be  remembered  that 
there  is  no  unique  definition  of  the  number  of  electrons  that  are  associated  with  a 
given  atom  or  nucleus  in  a  molecule.  In  the  case  of  CH4,  the  Mulliken  population 
analysis  results  in  a  net  charge  on  hydrogen  which  is  double  that  given  by  the 
Lowdin  method  at  the  ST0-3G  level  (  0.06  and  0.03,  respectively  ).  The  basis  set 
used  also  has  a  profound  effect  on  the  net  charge  which  results.  A  Mulliken 
population  analysis  on  CH4  gives  hydrogen  a  net  charge  of  0.06  at  the  STO-3G 
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level  but  0.12  at  the  6-3 IG**  level.l  Despite  its  drawbacks,  the  Mulliken 
population  analysis  is  readily  obtained  and  has  been  shown  to  give  exceptional 
results. 

As  discussed  previously,  the  method  employed  by  Saeb0,  et  al.  was  simply 
to  optimize  the  free  molecule  surrounded  by  the  point  charges.  They  chose  to  place 
the  point  charges  in  positions  as  determined  by  X-ray  diffraction  and  to  assign  them 
net  charges  based  on  a  Mulliken  population  analysis  of  the  optimized  free  molecule. 
Despite  their  convincing  results,  they  suggested  two  improvements  within  the 
frame  of  their  model.  The  first  was  to  introduce  the  new  net  charges  calculated  for 
the  free  molecule  in  the  surrounding  point  charges.  Secondly,  they  recommended 
increasing  the  size  of  the  basis  seL^  These  questions  and  several  others  are  the 
focus  of  my  study. 

Based  upon  the  suggestions  of  Saeb0  and  coworkers,  the  expansion  of  the 
basis  set  by  inclusion  of  polarization  functions  is  a  logical  extension  of  their  work. 
The  increased  flexibility  offered  by  the  polarization  functions  is  an  improvement  in 
the  basis  set  over  that  used  in  their  study.  Also,  the  idea  of  introducing  the  new 
point  charge  values  calculated  from  the  free  molecule  suggests  an  iterative  approach 
to  me.  Although  Saeb0  and  coworkers  envisioned  only  one  additional  optimization 
cycle,  I  felt  that  multiple  optimization  cycles  would  allow  more  flexibility  in  the 
model  by  allowing  it  to  dynamically  change.  The  iterative  technique  was  applied  in 
my  study  by  optimizing  the  free  molecule  inside  the  point  charges  based  on  both  X- 
ray  and  neutron  diffraction  structures.  Then  the  net  charges  were  recalculated  by 
running  a  Mulliken  population  analysis  on  the  new  "free"  molecule  geometry.  The 
new  net  charges  were  assigned  to  the  experimental  point  charge  locations  and  the 
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cycle  was  repeated  until  the  free  molecule  converged  in  both  structure  and  energy. 
This  should  eUminate  some  of  the  error  associated  with  the  Mulliken  population 
analysis. 

One  problem  with  the  Saeb0  method  is  that  it  requires  a  detailed 
experimental  structure  of  the  soUd.  While  X-ray  diffraction  structures  are  available 
for  a  wide  range  of  compounds,  they  quite  often  are  unable  to  locate  the  hydrogens 
which  was  the  case  for  both  cyanoformamide  and  sulfuric  acid.  Then,  in  order  to 
implement  their  method,  the  positions  of  the  hydrogens  must  be  assumed  which 
certainly  could  have  profound  effects  on  the  results.  In  addition,  experimental 
results  are  prone  to  mistakes  just  like  any  other  method.  For  sulfuric  acid,  the  early 
X-ray  diffraction  results  determined  the  symmetry  to  be  C2V  instead  of  the  actual 
C2.  Therefore,  it  would  be  helpful  to  have  a  method  which  was  not  quite  so 
dependent  on  the  experimental  results. 

One  solution  is  to  position  the  point  charges  based  upon  the  free  molecule 
instead  of  the  experimental  structures.  The  coordinates  of  the  free  molecule  can 
easily  be  converted  from  Cartesian  to  fractional  coordinates  and  then  used  to 
generate  the  point  charges.  That  requires  only  the  cell  dimensions  and  the  space 
group  of  the  crystal  and  it  preserves  the  internal  geometry  of  the  free  molecule. 
These  parameters  can  be  determined  from  experiment,  assumed  from  similar 
crystals,  or  simply  guessed.  This  would  permit  calculations  on  molecules  for 
which  there  are  no  experimental  results  or  for  which  incomplete  crystal  structures 
exist.  However,  the  major  advantage  of  my  method  would  be  for  calculations  on 
molecules  where  the  electron  density  is  shifted  away  from  the  nuclear  center  and  for 
which  only  X-ray  diffraction  crystal  structures  exist.  One  of  the  disadvantages  of 
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X-ray  crystallography  is  the  fact  that  it  gives  bond  lengths  which  are  not  very 
accurate.  That  is  due  to  the  manner  in  which  X-rays  interact  with  matter.  The 
interatomic  distances  determined  by  X-ray  diffraction  represent  the  distance 
between  centers  of  electron  density  and  not  nuclei.  In  molecules  of  this  type.  X-ray 
diffraction  can  give  results  for  bond  lengths  which  greatly  differ  from  the  true  bond 
lengths.^  The  error  in  bond  length  could  cause  profound  effects  on  the  structure 
calculated  in  the  point  charge  model.  The  technique  of  basing  the  point  charges  on 
the  free  molecule  has  been  applied  in  my  study  using  both  the  X-ray  and  neutron 
diffraction  data.  In  these  calculations,  the  only  experimental  information  that  was 
used  were  the  cell  dimensions  and  the  space  group.  The  calculations  were  repeated 
in  an  iterative  manner  until  the  "free"  molecule  converged  in  both  structure  and 
energy.  By  basing  the  point  charges  on  the  "free"  molecules,  any  errors  caused  by 
inaccurate  X-ray  bond  lengths  should  be  eliminated. 

Since  the  value  of  net  charge  determined  by  the  Mulliken  population 
analysis  seems  to  be  somewhat  arbitrary,  the  structural  results  should  have  a 
dependence  on  the  charges  that  are  generated.  However,  it  is  not  known  whether 
this  is  a  significant  factor  affecting  the  calculated  structure.  The  dependence  on  the 
magnitude  of  the  charges  was  approached  in  two  manners.  First,  a  Mulliken 
population  analysis  was  run  on  both  the  experimental  X-ray  and  neutron  diffraction 
structures.  That  required  a  single  point  SCF  calculation  be  done  at  the  experimental 
geometry  in  order  to  generate  a  wave  function.  The  resulting  point  charges  varied 
by  as  much  as  15%  from  those  calculated  for  the  optimized  free  molecule.  Second, 
the  value  of  the  point  charges  was  normalized  to  one  of  the  previously  calculated 
charges.  This  preserved  the  neutrality  of  the  molecule  and  resulted  in  charges  for 
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the  atoms  about  double  those  for  the  free  molecule.  The  new  charges  were  than 
used  to  compute  geometries  which  were  compared  to  those  calculated  using  the 
charges  based  on  the  free  molecule. 

Although  the  increased  computation  time  required  by  the  addition  of  point 
charges  is  small,  on  the  order  of  500-700  seconds  per  run,  it  would  still  be 
beneficial  to  limit  the  point  charges  to  as  few  as  possible.  The  increased 
computation  time  becomes  a  significant  factor  for  the  iterative  methods  which 
require  up  to  30  runs.  In  the  Saeb0  study,  the  closest  neighbor  molecules  were 
included.  But  there  are  obviously  pieces  of  these  molecules  which  are  far  enough 
away  from  the  free  molecule  such  that  their  electrostatic  interaction  is  essentially 
zero.  Hence,  if  they  were  not  included,  the  computation  time  could  be  further 
reduced.  In  order  to  investigate  the  effect  on  the  computed  structure,  calculations 
were  made  where  only  atoms  within  a  given  sphere  from  the  center  were  included. 
Because  of  the  requirement  to  maintain  a  neutral  charge,  a  dummy  charge  was 
generated  to  make  the  sum  of  the  point  charges  zero  and  it  was  placed  in  excess  of 
150  A  from  the  free  molecule. 

CONVERGENCE  CRITERIA 

The  point  charge  calculations  require  a  different  type  of  convergence  criteria 
than  normal  calculations  on  isolated  molecules.  This  is  due  to  the  residual 
Cartesian  forces  which  result  from  the  presence  of  the  point  charges.  In  the  point 
charge  calculations,  the  optimized  geometry  occurs  when  the  internal  forces  are 
small  but  the  Cartesian  forces  may  still  be  quite  large. 

An  isolated  molecule  geometry  is  considered  optimized  when  the  following 
conditions  are  met:^^ 
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1 .  All  Cartesian  forces  are  less  than  0.005  mdyne/angstrom. 

2.  Stretching  and  bending  internal  coordinate  forces  are  less  than 
0.005  mdyne/angstrom  or  /radian. 

3.  Torsional  internal  coordinate  forces  are  less  than  0.001 
mdyne/radian. 

4.  The  internal  force  on  torsional  coordinates  has  changed  sign. 
Conditions  1  through  3  are  the  normal  conditions  used  to  determine 

geometry  convergence.  Condition  4  is  absolutely  essential  to  ensure  convergence 
of  low  energy  torsional  and  out-of-plane  motions.34  The  above  conditions  should 
give  convergence  of  calculated  bond  lengths  to  ±  0.002  A,  bond  angles  to  ±  0.2°, 
and  torsional  angles  to  ±  2.0°. 

As  discussed  earlier,  the  convergence  criteria  for  the  point  charge 
calculations  must  be  slightly  different.  For  these  calculations,  conditions  2  through 
4  above  must  be  met  with  the  additional  requirement  for  all  Cartesian  coordinate 
displacements  to  be  less  than  0.0005  A.  The  restriction  on  the  change  in  Cartesian 
coordinates  is  based  upon  the  work  by  Saeb0  and  my  review  of  other  isolated 
molecule  optimizations  when  conditions  1  through  4  were  met. 


CHAPTER  THREE 

SULFURIC  ACID 

One  reason  for  the  great  interest  in  sulfuric  acid  is  because  of  its  presence  in 
"acid"  rain.  The  harmful  effects  of  "acid"  rain  have  been  widely  reported  over  the 
last  decade.  These  detrimental  effects  include  large  fish  kills,  contamination  of 
water  supplies,  and  damage  to  automobile  finishes.  The  Northeastern  portion  of 
the  United  States  has  been  hit  especially  hard.  Research  has  found  this  to  be 
predominantly  due  to  the  increased  burning  of  high  sulfur  coal  in  the  Midwest. 
More  recentiy,  the  effects  of  "acid"  rain  have  been  reported  worldwide.  Many  old 
churches  and  monuments  in  Europe  are  rapidly  being  destroyed  by  "acid"  rain.  For 
these  reasons  and  several  others,  such  as  the  detection  of  sulfuric  acid  in  the 
Venusian  atmosphere,  many  studies  have  been  done  on  sulfuric  acid.  These 
studies  have  centered  on  the  formation  and  structure  of  sulfuric  acid,  especially  the 
intermediates  that  are  involved.  Results  from  some  of  these  studies  have  lead  to  an 
understanding  of  the  formation  of  sulfuric  acid  and  its  detection  in  the  vapor  phase. 

As  discussed  in  Chapter  One,  sulfuric  acid  undergoes  significant  structural 
changes  transforming  from  a  gas  to  a  solid.  Regardless  of  the  phase,  the  symmetry 
remains  C2.  The  H2SO4  molecule  with  C2  symmetry  is  shown  in  Figure  3-1  along 
with  the  atom  designations  that  will  be  used  throughout  the  paper.  These  structural 
alterations  can  easily  be  seen  in  a  comparison  of  the  gas  and  solid  structures  as 
shown  in  Table  3-1.  The  most  dramatic  change  is  the  shortening  of  the  S-0  single 
bond  on  the  order  of  0.05  A.  In  addition,  there  is  a  noticeable  decrease  in  the 
0=S=0  bond  angle  coupled  with  an  increase  in  the  O-S-O  bond  angle.  These 
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Figure  3-1 

Atom  Numbering  Scheme  for  Sulfuric  Acid  (H2SO4) 
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Table  3-la 

Sulfuric  Acid  (H2SO4)  Gas  and  Crystal  Geometries 

MWb        X-RAY(1978)c         NDd         X-RAY(1965)e 

solid 

1.528(5) 

1.419(5) 

f 

f 

103.7(4) 

118.3(4) 

110.6(3) 

106.3(3) 

f 

f 

a.  Bond  lengths  in  angstroms.  Bond  angles  in  degrees. 

b.  Preferred  microwave  structural  parameters  and  uncertainties  from  reference  35. 

c.  X-ray  diffraction  geometry  with  estimated  standard  deviations  (esd's)  in 

parentheses  from  reference  36. 

d.  Neutron-diffraction  10°K  structure  from  reference  37.  The  estimated 

uncertainties  are  reported  in  the  paper  as  ±0.05-0.10  A. 

e.  X-ray  diffraction  geometry  with  esd's  in  parentheses  from  reference  38. 

f.  The  hydrogens  were  not  located  in  either  X-ray  diffraction  study. 

g.  T(ijkl)  is  the  torsional  angle  between  the  ijk  and  jkl  planes. 


FHase 

gas 

S-Ol 

1.57410.01 

S=02 

1.422±0.01 

Ol-H 

0.97±0.01 

ZH-Ol-S 

108.5±1.5 

zoi-s-or 

101.3±1 

Z02-S-02' 

123.3±1 

ZOl-S-02 

108.6±0.5 

Z01-S-02' 

106.410.5 

T(H-01-S-01')g 

90.9±1 

x(H-01-S-02)g 

-20.8±1 

solid 

solid 

1.48 

1.535(15) 

1.49 

1.426(15) 

0.91 

f 

113.6 

f 

112.8 

104.0(10) 

110.1 

118.6(10) 

108.0 

110.5(10) 

109.0 

105.9(10) 

84.1 

f 

-36.5 

f 
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significant  changes  are  attributed  to  the  hydrogen  bond  formation  in  the  crystal.^  I 
have  been  able  to  reproduce  these  structural  changes  with  the  point  charge  model. 

The  remainder  of  this  chapter  will  deal  with  the  results  I  have  obtained  for 
both  the  gas  and  solid  phases  of  sulfuric  acid.  The  gas  phase  will  be  discussed 
first,  followed  by  the  crystal  phase. 

GAS 

Although  vapor  pressure  studies  in  1923  showed  that  solutions  with 
concentrations  exceeding  85%  had  measurable  vapor  pressures  above  200°  F,  early 
studies  of  gaseous  sulfuric  acid  were  hampered  by  several  problems.  An  infrared 
(IR)  spectrum  of  gaseous  sulfuric  acid  was  not  completed  until  1965,  because  of  its 
low  volatility  and  attack  on  windows  at  high  temperatures.  In  that  study,  the 
observed  H2SO4  bands  were  assigned  in  agreement  with  the  expected  chemical 
series  (X-SO2-Y:  X,Y=  F,  OH,  CI, . .  .,CH3).  Based  on  those  assignments,  the 
researchers  incorrecdy  concluded  that  free  H2SO4  molecules  had  an  approximate 

tetrahedral  configuration  with  C2V  symmetry  similar  to  the  other  X-SO2-Y 
compounds.39 

The  first  molecular  orbital  calculation  of  the  gas  phase  structure  of  H2SO4 
was  undertaken  in  1978.  The  purpose  of  that  study  was  to  gain  insight  into  the  gas 
phase  reaction  of  SO3  and  H2O  to  form  H2SO4  using  semi-empirical  CNDO/2 
calculations.  Those  researchers  incorrecdy  calculated  the  symmetry  of  H2SO4  as 
C2V  in  agreement  with  the  earlier  IR  study.  The  calculated  structural  parameters 
were  01-H=  1.04  A,  S-01=  1.63  A,  S-02=  1.56  A,  Z01-S-01'=  97°,  ZH-01-S= 
97.5°,  and  Z02-S-02*=  116°  witii  the  H-01-S-01*-H'  atoms  being  planar.^O 
Their  configuration  is  referred  to  as  the  compact  C2V  geometry  where  the  OH 
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groups  point  toward  and  bisect  the  terminal  SO2  group  with  a  torsional  angle  (ZH- 
Ol-S-Or  and  ZH'-Ol'-S-Ol)  of  180°. 

The  first  ab  initio  calculation  on  H2SO4  was  completed  in  1980  due  to  the 
importance  of  sulfur-oxygen  compounds  in  air  pollution  chemistry,  specifically 
"acid"  rain.^l  The  optimum  geometries  in  that  study  were  calculated  using  the 
ST0-3G*  basis  set.  Previous  calculations  on  similar  compounds  using  the  STO- 
3G  or  4-3 IG  basis  sets  had  shown  poor  agreement  between  experiment  and 
calculated  structural  parameters  though  the  STO-3G*  basis  set  had  been  shown  to 
give  acceptable  geometries.  Since  the  researcher's  principal  interest  was  in  the 
geometry  about  the  sulfur  atom,  the  OH  length  and  SOH  angle  were  fixed  at  0.994 
A  and  105°,  respectively  [these  values  corresponded  to  their  optimum  values  in 
S02(OH)].  The  calculated  STO-3G*  geometry  for  sulfuric  acid  was  S-01=  1.63 
A,  S-02=1.45  A,  Z01-S-01'=  99°,  Z01-S-02=  107°,  and  Z02-S-02'=  125°. 
The  H-O-S-O  dihedral  angles  were  set  to  90°  by  the  researchers  based  upon  the  best 
ideahzed  conformation.  The  addition  of  d  functions  to  the  basis  set  was  shown  to 
greatly  improve  the  calculated  energy.  The  energies  calculated  at  the  optimum 
geometry  for  the  ST0-3G*  and  4-3 IG*  basis  sets  showed  an  energy  improvement 
of  316.1  and  170.6  kcal/mole,  respectively,  over  the  energies  at  the  STO-3G  and  4- 
31Glevels.4l 

The  correct  C2  symmetry  was  finally  determined  in  1980  from  the 
microwave  spectra  of  four  isotopic  species  of  H2SO4.  From  those  spectra,  the 
researchers  determined  a  detailed  structure  for  gaseous  H2SO4  which  is  shown  in 
Table  3-1.  In  contrast  to  the  previously  determined  compact  C2V  geometry,  the  OH 
bonds  were  found  to  have  rotated  past  the  S02  and  S02'  bonds  resulting  in 
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dihedral  angles  of  -20.8°  with  them.  That  structure  possessed  a  H-Ol-S-Ol' 
dihedral  angle  of  90.9°  as  compared  to  180°  found  in  the  CNDO/2  calculation.35 
The  H2SO4  molecule  with  C2  symmetry  is  shown  in  Figure  3-1. 

The  continued  interest  in  H2SO4  prompted  two  additional  ab  initio  studies 
of  its  structure  and  torsional  modes.'*2,43  Both  of  these  studies  determined  the 
potential  energy  surface  for  H2SO4  as  a  function  of  <^  (ZH-Ol-S-OT)  and  (()'  (Z 
H'-Ol'-S-Ol).  The  p,p  structure  ("p"  denotes  periplanar  or  cis)  is  shown  in 
Figure  3-2  where  (()=(1)'=0°.  These  researchers  performed  most  of  their  calculations 
at  the  STO-3G  level  except  for  some  selected  conformations  which  were  calculated 
at  the  4-3 IG  level.  Despite  the  different  structures  used  for  sulfuric  acid,  as 
discussed  below,  the  results  of  these  two  studies  were  surprisingly  similar. 

The  first  study  was  started  prior  to  the  completion  of  the  microwave  study. 
As  a  result,  Kaliannan  and  coworkers  used  the  geometry  of  dimethyl  sulfate 
[S02(OCH3)]  obtained  by  electron  diffraction  in  their  study  of  sulfuric  acid.  The 
0-H  bond  length  and  H-O-S  bond  angle  were  set  at  0.96  A  and  109.47°, 
respectively.  Their  construction  of  the  potential  energy  surface  consisted  of 
calculations  at  22  conformations  of  <^  and  (|)'.  They  found  the  global  minimum 
occurring  at  (})=(j)'=90°  in  good  agreement  with  the  microwave  study  value  of 
({)=(j)'=90.9°.  Their  calculated  dipole  moment  at  that  conformation  was  2.51  debye 
(D)  compared  to  the  microwave  value  of  2.725  D.  A  secondary  minimum  was 
found  to  occur  at  (})=-(()'= 120°  with  an  energy  about  1.2  kcal/mole  above  the  global 
minimum.  They  also  determined  that  the  variation  in  energy  was  very  small  in  the 
range  90°-120°  for  both  (|)  and  (])',  suggesting  that  the  molecule  is  reasonably  flexible 
in  that  range.  A  comparison  of  their  potential  energy  surface  to  that  of  H2PO4" 
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Figure  3-2 
H2SO4  p,p  structure  ({)=(()'=0°  ("p"  denotes  periplanar  or  cis) 


or 
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found  the  two  to  be  very  similar  with  the  global  minimum  occurring  at  (t)=(t)'=90°  in 
both  cases.  The  results  of  their  4-3 IG  calculations  at  the  STO-3G  optimized 
geometries  confmned  the  presence  of  the  global  minimum  at  ^=(^'=90°  in  agreement 
with  the  microwave  experiment.^^ 

Lohr  has  also  conducted  a  study  of  the  structure  and  torsional  modes  of 
H2SO4.  In  his  study,  the  microwave  values  for  bond  lengths  and  bond  angles 
were  used.  He  constructed  his  potential  surface  from  43  calculated  points  at  the 
ST0-3G  level.  Lohr  found  the  global  minimum  had  C2  symmetry  and  occurred  at 
({)=(j)'=98.6°,  somewhat  larger  than  the  microwave  value  of  <j)=(t) -90.9°  .  A  local 
minimum  with  Cs  symmetry  was  also  found  at  <])=-(|)'=104°  with  an  energy  of  1.36 
kcal/mole  relative  to  the  global  minimum.  Calculations  at  eleven  selected  points 
were  done  using  the  larger  4-3 IG  basis  set.  The  structure  with  C2  symmetry  was 
found  to  have  a  minimum  energy  at  (j)=(|)'=94.9°  which  is  somewhat  closer  to  the 
microwave  value  than  is  the  STO-3G  value.'^^  The  results  of  both  torsional  mode 
studies  confirm  the  C2  symmetry  which  is  observed  in  the  microwave  study. 

I  optimized  the  gas  phase  of  sulfuric  acid  with  several  different  basis  sets. 
The  programs  TEXAS,  TEXAS(D),  and  GAUSSIAN-82  44  were  used  to  obtain  an 
optimum  structure  at  the  6-3 IG**  level.  In  addition,  calculations  were  conducted 
at  the  4-3  IG**  level  using  TEXAS  and  at  the  STO-3G*  level  using  GAUSSIAN- 
85  45.  The  results  of  my  calculations  on  gaseous  sulfuric  acid  are  presented  in 
Table  3-2. 

The  6-3 IG**  basis  set  admirably  reproduces  the  experimental  structure. 
With  very  few  exceptions,  the  calculated  6-3 IG**  structure  is  within  the 
experimental  uncertainties  of  the  microwave  study.  The  major  variances  occur  in 
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Table  3-2^ 
Sulfuric  Acid  (H2SO4)  Gas  Phase  Geometries 


MW^ 


5.31G**c      6-3 IG**'^      6-31G**<^      4-310**'^     STO-3G*^ 


TEXAS(D) 

G-82 

G-85 

S-01 

1.574±0.01 

1.570 

1.569 

1.569 

1.563 

1.620 

S=02 

1.422±0.01 

1.411 

1.411 

1.411 

1.408 

1.446 

01-H 

0.97±0.01 

0.950 

0.950 

0.950 

0.950 

0.993 

ZH-Ol-S 

108.511.5 

110.7 

110.8 

110.8 

111.2 

105.6 

ZOl-S-01' 

101.3±1 

101.8 

101.8 

101.8 

101.7 

100.1 

Z02-S-02' 

123.3±1 

123.5 

123.6 

123.6 

123.5 

125.6 

ZOl-S-02 

108.6±0.5 

108.1 

108.1 

108.1 

108.3 

108.9 

ZOl-S-02' 

106.4±0.5 

106.7 

106.7 

106.6 

106.5 

105.2 

tCH-OI-S-OI") 

90.9±1 

88.9 

89.2 

88.8 

87.8 

79.1 

T(H-01-S-02) 

-20.8±1 

-23.2 

-22.9 

-23.3 

-24.2 

-31.0 

Energy: 

Haitree  -698.051422  -698.052804  -689.766048 

-698.052805  -697.659157 

kcal/mole  -438,033.781  -438,034.65  -432,834.63 

-438,034.649  -437,787.631 

a.  Bond  lengths  in  angstroms.  Bond  angles  in  degrees. 

b.  Preferred  microwave  structural  parameters  and  uncertainties  from  reference  35. 

c.  This  work.  Calculations  run  with  TEXAS  unless  noted  otherwise. 
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the  torsional  angles  which  are  notoriously  difficult  to  reproduce.  Considering  the  ± 
2.0°  uncertainty  in  the  value  of  the  calculated  torsion  angles,  the  agreement  is  very 
good. 

The  development  of  TEXAS  (D)  by  Pulay's  group  prompted  me  to  optimize 
the  gas  phase  structure  with  it  for  comparison.  As  shown  in  Table  3-2,  the  6- 
31G**  structures  are  identical  for  both  TEXAS  and  TEXAS(D).  My  calculations 
confirm  that  the  displaced  p  functions  used  to  construct  d  functions  in  TEXAS  are 
more  than  adequate.  The  only  difference  between  these  two  calculated  structures  is 
in  the  energy.  The  TEXAS(D)  energy  is  0.868  kcal/mole  lower  than  that  calculated 
in  TEXAS.  The  improved  energy  in  TEXAS (D)  must  be  due  to  the  use  of  true  d 
functions.  For  most  calculations,  this  energy  difference  is  insignificant.  Based  on 
calculated  structures  and  energy,  neither  TEXAS  nor  TEXAS(D)  offer  any 
advantage  over  the  other  with  the  exception  of  the  slightiy  lower  energy  in 
TEXAS (D).  In  terms  of  computation  time,  TEXAS (D)  is  superior,  posting 
decreases  of  up  to  30%  over  TEXAS.  TEXAS (D),  however,  requires  double  the 
amount  of  disk  space  to  store  the  integrals  as  compared  to  TEXAS. 

Because  of  the  recent  acquisition  of  a  Cray  version  of  GAUSSIAN-82, 1 
decided  to  make  a  comparison  between  it  and  TEXAS (D).  Since  both  programs 
incorporate  true  d  functions,  I  expected  identical  results.  As  shown  in  Table  3-2, 
the  results  are  essentially  identical.  The  few  tenths  of  a  degree  difference  between 
torsional  angles  is  well  within  the  calculated  uncertainty  of  ±  2.0°.  The  TEXAS(D) 
program  calculated  an  energy  which  is  only  6x10"^  kcal/mole  lower  than 
GAUSSIAN-82.  Based  upon  my  comparison,  both  TEXAS(D)  and  GAUSSIAN- 
82  do  give  identical  results  for  calculations  conducted  with  the  same  basis  set. 
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Initially,  I  was  interested  in  reducing  the  computation  time  required  for  one 
optimization  run.  As  a  result,  I  completed  an  optimization  run  of  gaseous  sulfuric 
acid  at  the  4-31G**  level.  The  4-31G**  basis  set  I  utilized  contained  a  4-31G** 
basis  on  sulfur  and  6-3 IG**  on  both  oxygen  and  hydrogen.  As  shown  in  Table  3- 
2,  the  4-3 IG**  basis  set  gave  structural  parameters  almost  identical  to  those 
computed  with  TEXAS  at  the  6-3 IG**  level.  The  major  distinctions  are  the  shorter 
Ol-S  bond  length  and  the  torsional  angles  at  the  4-3  IG**  level.  The  computation 
time  saved  at  the  4-3  IG**  level  was  significant.  The  4-3 IG**  calculation  ran  15% 
faster  than  the  6-31G**  calculation,  1761  seconds  compared  to  2066  seconds.  The 
time  savings  is  due  to  the  reduction  in  the  number  of  primitive  Gaussian  functions. 
Although  the  number  of  contracted  Gaussians  is  the  same  for  both  basis  sets  at  89, 
the  number  of  primitive  Gaussians  is  reduced  from  208  to  198  in  going  from  the  6- 
31G**  basis  set  to  the  4-31G**  basis  set.  The  major  disagreement  between  the 
two  calculations  is  the  total  energy.  The  6-3 IG**  calculation  results  in  an  energy 
247  kcal/mole  lower  than  the  energy  at  the  4-3  IG**  level.  For  many  calculations, 
the  4-31G**  basis  set  would  be  satisfactory.  Since  I  was  concerned  with 
observing  structural  changes,  I  chose  to  use  the  more  accurate  6-3 IG**  basis  set  in 
the  crystal  calculations.  In  addition,  the  6-3 IG**  basis  set  ((lls4pld/4slp) 
contracted  to  [4s2pld/2slp] }  represents  a  significant  improvement  over  the  basis 
set  used  by  Saeb0  and  coworkers  {(7s3p/4s)  [4s2p/2s]}.  The  6-3 IG**  basis  set 
not  only  increases  the  number  of  primitive  functions  but  also  adds  polarization 
functions  to  hydrogen  and  tiie  heavy  atoms. 

In  order  to  make  comparisons  with  some  dimer  calculations,  I  optimized 
sulfuric  acid  at  the  ST0-3G*  level.  As  shown  in  Table  3-2,  the  calculated 
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parameters  at  this  level  differ  greatly  from  the  experimental  structure.  The  large 
error  in  the  torsional  angles  at  the  STO-3G*  level  confirms  the  need  for  polarization 
functions  on  oxygen  to  accurately  calculate  these  angles.  More  importantly,  while 
no  symmetry  restrictions  were  placed  on  any  of  the  other  calculations,  I  had  to 
restrict  the  STO-3G*  calculation  to  C2  in  order  to  prevent  the  program  from 
determining  a  structure  with  Ci  symmetry.  While  minimal  basis  set  calculations 
can  yield  valuable  information,  the  need  for  accurate  structural  information  almost 
always  requires  using  a  larger  basis  set. 

CRYSTAL 

The  first  X-ray  diffraction  study  of  sulfuric  acid  was  completed  in  1954  by 
Pascard.  That  study  incorrecdy  concluded  that  sulfuric  acid  crystallized  in  the 
noncentrosymmetric  space  group  C2  and  possessed  C2V  symmetry .^^6  jn  1964, 
Pascard-Billy  refined  the  original  intensity  data  with  anisotropic  temperature  factors 
for  oxygen  atoms  to  a  reliability  (R)  factor  of  only  0.124.  The  reported  dimensions 
of  the  H2SO4  molecule  deviated  significantly  from  idealized  C2V  symmetry. 
Pascard-Billy  determined  that  the  correct  space  group  was  actually  C2/c  and  the 
H2SO4  molecule  possessed  C2  symmetry.  Although  Pascard  was  not  able  to  locate 
the  hydrogens,  he  theorized  that  they  each  formed  a  hydrogen  bond  with  a  double 
bonded  oxygen  on  an  adjacent  molecule.  In  that  configuration,  one  H2SO4 
molecule  is  hydrogen  bonded  to  four  adjacent  molecules  with  the  hydrogen  atom 
lying  on  the  -01 . .  .02=  hne.38 

In  1978,  Yu  and  Mak  felt  that  the  atomic  parameters  reported  earlier  could 
be  further  refined  by  full-matrix  anisotropic  least  squares  with  correction  for 
secondary  extinction.  Their  study  showed  significant  improvement  of  the  standard 
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deviations  in  atomic  parameters  over  those  reported  earlier  with  an  R  value  of 
0.099.  Their  measured  hydrogen  bond  length  (Ol-H. .  .02)  of  2.624  A  belongs  to 
the  short  type.  They  concluded  that  the  hydrogen  atom  almost  lay  on  the 
intemuclear  Ol . .  .02  line  if  a  tetrahedral  S-O-H  angle  was  assumed.36 

The  hydrogen  atoms  were  finally  located  in  1983  by  Moodenbaugh  and 
coworkers  using  powder  neutron-diffraction  techniques.  Diffraction  data  was 
collected  near  the  melting  point  (240°K)  and  at  10°K  to  search  for  order-disorder 
transitions  in  the  hydrogen  bond  network.  The  neutron-diffraction  scans  revealed, 
to  the  limits  of  detection,  that  the  sulfuric  acid  sample  was  single  phase  with 
monoclinic  symmetry  at  both  temperatures.  The  principal  conclusions  from  their 
study  were  that  the  hydrogens  were  ordered  at  both  temperatures  and  the 
parameters  of  the  other  atoms  were  qualitatively  consistent  with  the  single-crystal 
X-ray  studies.^^ 

The  generation  of  the  point  charges  used  in  my  calculations  was  faciUtated 
by  the  use  of  several  computer  programs.  The  BMW  1^7  (Best  Molecular  Fit) 
program  derives  the  best  least-squares  fit  between  two  sets  of  atoms.  I  used 
BMFIT  in  order  to  convert  calculated  structures  into  fractional  coordinates  in  the 
crystal  space.  The  program  DAESD^S  calculates  interatomic  distances  and  angles 
in  crystal  structures.  The  DAESD  program  was  used  to  generate  all  the  molecules 
within  a  given  unit  cell.  The  programs  ZERO  and  MOLE  were  written  by  me  in 
order  to  generate  the  actual  pomt  charge  locations.  Both  programs  convert  the 
contents  of  a  unit  cell  into  a  molecule  and  its  nearest  neighbors.  In  the  program 
ZERO,  all  neighbor  atoms  within  a  given  sphere  of  radius  r  from  the  central 
molecule  are  punched.  In  addition,  a  dummy  charge  is  generated  and  placed  at 
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(100,100,100)  in  order  to  maintain  the  neutrality  of  the  point  charges.  In  MOLE, 
aU  neighbor  molecules  that  have  their  central  atom  within  r  from  the  central 
molecule  are  punched.  In  the  case  of  MOLE,  the  neutrality  of  charge  is  maintained 
by  including  whole  molecules  instead  of  individual  atoms. 

My  studies  have  shown  that  great  care  must  be  taken  in  generating  the  point 
charges.  As  discussed  in  Chapter  Two,  I  initially  tried  to  limit  the  number  of  point 
charges  in  order  to  reduce  the  computation  time.  I  did  this  by  using  the  program 
ZERO  to  generate  the  point  charges.  The  program  ZERO  allowed  me  to  include 
only  the  nearest  atoms  and  not  the  nearest  molecules.  I  could  then  eliminate  atoms 
on  adjacent  molecules  that  were  far  enough  away  that  they  had  essentially  zero 
electrostatic  interactions  with  the  central  molecule.  However,  this  method  proved 
to  be  entirely  unsatisfactory.  Not  only  was  the  savings  in  computation  time 
negligible,  but  this  method  also  introduced  artificial  symmetries  into  the  calculated 
structure.  In  my  particular  case,  the  point  charges  generated  using  the  program 
ZERO  artificially  reduced  the  symmetry  of  H2SO4  from  C2  to  Ci.  As  a  result,  all 
of  my  crystal  calculations  were  made  using  the  nearest  neighbor  molecules 
generated  by  MOLE.  After  careful  scrutiny  of  the  point  charges,  I  found  that 
including  the  18  nearest  H2SO4  molecules  insured  that  all  nearest  neighbor  atoms 
were  included  in  the  point  charges.  All  of  the  reported  calculated  crystal  structures 
used  126  point  charges  based  on  the  18  nearest  H2SO4  molecules. 

The  crystal  calculations  which  I  have  completed  can  be  divided  into  four 
categories  -  SF,  SE,  MP,  and  IF.  These  designations  are  used  throughout  the 
remainder  of  the  paper  and  are  defined  as  follows: 
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SF-       A  single  optimization  cycle  was  run.  (Note:  An  optimization  cycle 
means  that  the  convergence  criteria  of  Chapter  Two  have  been  met.) 
The  point  charge  magnitudes  are  based  on  the  optimized  free  (gas 
phase)  molecule.  The  point  charge  locations  are  based  on  the 
experimental  crystal  structures.  This  is  the  method  used  by  Saeb0 
and  coworkers. 

SE-       This  is  identical  to  the  SF  method  except  that  the  point  charge 

magnitudes  are  based  on  a  calculation  at  the  experimental  crystal 
geometry. 

MF-      This  method  is  a  continuation  of  the  SF  method  involving  multiple 
optimization  cycles.  In  this  case,  a  Mulliken  population  analysis  is 
run  on  the  new  "free"  molecule  geometry  obtained  from  the  SF 
method.  (Note:  The  term  "free",  in  quotation  marks,  refers  to  the 
central  molecule  surrounded  by  the  point  charges.)  The  new  net 
charges  are  assigned  to  the  experimental  point  charge  locations  and 
another  optimization  cycle  is  completed  on  the  "free"  molecule.  The 
process  is  repeated  until  the  "free"  molecule  converges  in  both 
structure  and  energy. 

IF-        This  is  an  iterative  method  involving  multiple  optimization  cycles 
similar  to  the  MF  method.  In  this  case,  however,  both  the  point 
charge  magnitudes  and  locations  are  based  on  the  optimized  free 
(gas  phase)  molecule.  The  only  information  used  from  the 
experimental  crystal  study  is  the  unit  cell  dimensions  and  space 
group.  After  each  optimization  cycle,  the  magnitude  and  location  of 
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the  point  charges  are  adjusted  to  those  corresponding  to  the  newly 
optimized  "free"  molecule.  The  process  is  continued  until  the  "free" 
molecule  converges  in  both  structure  and  energy. 

The  SE  method  was  conceived  by  me  to  test  the  sensitivity  of  the  calculated 
crystal  structures  to  the  magnitude  of  the  point  charges.  A  Mulliken  population 
analysis  was  run  on  the  wave  function  from  a  single  point  SCF  calculation  at  the 
experimental  geometry  to  obtain  the  experimental  net  charges.  As  expected,  the  net 
charges  from  the  experimental  geometry  differed  significantiy  from  those  of  the 
optimized  free  molecule.  Although  using  the  net  charges  from  the  free  molecule 
seems  logical  since  a  wave  function  already  exists,  no  study  has  been  made 
previously  to  compare  the  effects  of  the  two  sets  of  net  charges.  Since  a  single 
point  SCF  calculation  is  relatively  small  (about  650  seconds),  the  calculation  of  the 
experimental  net  charges  and  their  subsequent  use  in  the  SE  method  seems  a 
reasonable  altemative  to  the  SF  method. 

The  use  of  the  iterative  method  in  the  point  charge  model  is  an  extension  of 
an  improvement  proposed  by  Saeb0  and  coworkers.  They  simply  suggested  that 
the  point  charge  model  might  be  improved  by  including  the  net  charges  of  the 
newly  optimized  "free"  molecule.  My  idea  was  to  iterate  this  process  until  the 
molecule  converged  in  both  structure  and  energy.  My  iterative  technique  is 
therefore  self-consistent.  The  major  advantage  of  my  method  is  that  it  adds 
flexibility  to  the  point  charge  model.  A  significant  concern  with  the  point  charge 
model  is  that  it  may  not  incorporate  all  the  forces  which  bind  crystals  together.  As 
discussed  in  Chapter  Two,  all  of  the  binding  forces  in  crystals  involve  charged 
particle  interactions.  Since  the  model  basically  allows  for  electrostatic  interactions 
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between  the  electrons  and  the  point  charges,  the  iterative  method  should  accurately 
recreate  the  environment  inside  the  crystal.  In  addition,  any  effect  of  the  arbitrary 
assignment  of  net  charge  in  the  population  analysis  should  be  greatly  minimized 
since  the  net  charges  are  varied  in  response  to  their  environment. 

I  devised  the  IF  method  for  a  couple  of  reasons.  First,  experimental  crystal 
data  is  often  incomplete  or  inaccurate.  As  discussed  previously,  X-ray 
crystallography  frequently  cannot  locate  hydrogen  atoms.  Even  when  all  the  atoms 
are  located,  the  precision  with  which  an  atom  can  be  positioned  in  a  crystal  is 
subject  to  systematic  errors.  These  errors  can  be  divided  into  two  classes:  errors  in 
the  model  and  errors  in  the  data.  If  it  was  possible  to  collect  data  with  no  errors, 
many  significant  structural  errors  could  be  introduced  if  an  improper  model  is  used 
in  the  refinement.  The  IF  method  is  not  as  dependent  on  the  experimental  data 
since  it  uses  only  the  cell  dimensions  and  space  group  of  the  crystal.  The  cell 
dimensions  and  space  group  are  more  easily  determined  and  can  be  found  even  if 
all  the  atoms  are  not  located.  Second,  a  principal  concern  with  the  point  charge 
model  is  that  it  may  not  adequately  represent  the  crystal  environment.  The  IF 
method  allows  for  movement  of  the  point  charges.  Therefore,  the  point  charges 
should  ultimately  reside  in  the  crystal  locations  if  the  point  charges  properly  model 
the  crystal  environment.  Thus,  the  IF  method  serves  as  a  test  of  the  validity  of  the 
point  charge  model. 

The  first  calculations  that  I  completed  were  based  on  the  X-ray  study  of  Yu 
and  Mak.  Since  the  hydrogens  were  not  located  in  their  study,  I  assigned  their 
location  by  assuming  an  Ol-H  bond  length  of  0.95  A  and  a  hydrogen  bond  angle 
(ZOl-H. .  .02)of  180°.  The  results  of  these  calculations  are  presented  in  Table  3-3. 
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In  addition  to  the  calculations  based  on  the  X-ray  study,  I  also  conducted 
calculations  based  on  the  neutron-diffraction  study.  The  10°K  neutron-diffraction 
data  was  used  as  the  basis  for  the  calculations  because  of  its  minimal  thermal 
motion  and  due  to  errors  in  the  240°K  data  I  discovered.  The  results  of  these 
calculations  are  presented  in  Table  3-4.  The  numerous  problems  with  the  neutron- 
diffraction  study  will  be  discussed  later. 

One  aspect  that  all  of  the  calculated  solid  structures  share  is  a  lower  energy 
than  the  gas  phase.  The  energy  improvement  in  the  solid  structures  is  on  the  order 
of  50-70  kcal/mole.  The  lower  energy  in  the  crystal  is  expected  based  on  the 
stabilizing  effect  of  the  neighbor  molecules  in  the  solid.  The  calculated  energy 
improvement  which  corresponds  to  12.5-17.5  kcal/mole  per  hydrogen  bond 
compares  favorably  to  7.8  kcal/mole  for  cyanoformamide^,  5.4  and  5.9  kcal/mole 
for  CH3CN...HOCH3,  respectively^^^  3.9-9.6  kcal/mole  for  the  water  dimer 
(depending  upon  basis  set)^^,  and  7.0-10.3  kcal/mole  for  the  various 
monohydrated  formamide  complexes^l.  In  addition,  all  of  the  calculated  crystal 
structures  mimic  the  major  structural  changes  which  occur  in  H2SO4  in  going  from 
a  gas  to  a  solid.  Specifically,  the  calculated  crystal  structures  show  a  substantial 
shortening  of  the  S-0  bond  along  with  an  increase  in  the  ZOl-S-Ol'  and  a 
decrease  in  the  Z02-S-02'  compared  to  the  gas.  Further,  aU  of  the  calculated 
crystal  structures  agree  with  the  general  and  very  characteristic  structural  feature 
between  bond  angles  for  sulphones  (XSO2Y).  Hargittai  has  demonstrated  that  the 
ZX-S-Y  <  ZX(Y)-S=0  <  Z  0=S=0,  regardless  of  the  X  and  Y  ligands.  The 
Z.  0=S=0  is  considerably  larger  than  the  ideal  tetrahedral  angle  while  the  Z  X-S-Y 
is  considerably  smaller.52  While  tiiese  results  alone  substantiate  the  point  charge 
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Table  3-3^ 

Sulfuric  Acid  (H2SO4)  X-ray  Crystal  Geometry 


X 

3^Y(1978; 

)b     spc 

SEC 

Mb^ 

Ipc 

s-01 

1.528(5) 

1.558 

1.557 

1.553 

1.552 

S=02 

1.419(5) 

1.422 

1.424 

1.428 

1.428 

Ol-H 

0.950d 

0.971 

0.972 

0.978 

0.978 

01...02e 

2.624(8) 

2.553 

2.553 

2.550 

2.556 

ZH-Ol-S 

109.9d 

113.9 

114.1 

114.7 

114.8 

zoi-s-or 

103.7(4) 

104.3 

104.6 

105.3 

105.2 

Z02-S-02' 

118.3(4) 

120.8 

120.5 

119.6 

119.4 

ZOl-S-02 

110.6(3) 

107.3 

107.3 

107.2 

107.4 

Z01-S-02' 

106.3(3) 

108.0 

108.1 

108.4 

108.3 

ZS-01...02e 

110.0(3) 

111.3 

111.3 

111.3 

111.1 

Z01...02e-Se 

157.9(4) 

160.7 

160.7 

160.5 

160.4 

X(H-Ol-S-Ol') 

86.3d 

92.3 

92.0 

91.6 

89.5 

T(H-01-S-02) 

-27.3d 

-22.2 

-22.6 

-23.7 

-25.7 

Energy: 

Hartree 

-698.130620 

-698.155672 

-698.134681  -698.161332 

kcal/mole  -438,083.479  -438,099.199 

-438,086.027  -438,102.750 

a.  Bond  lengths  in  angstroms.  Bond  angles  in  degrees. 

b.  X-ray  diffraction  geometry  with  esd's  in  parentheses  from  reference  36. 
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c.  This  work.  Calculated  using  the  6-3 IG**  basis  set. 

d.  Based  on  setting  the  Ol-H  bond  length  to  0.95  A  and  the  Ol-H. .  .02  bond 
angle  to  180°. 

e.  Atom  is  located  at  fractional  coordinates  -l/2-i-x,l/2+y,z  on  an  adjacent 
molecule. 
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Table 

3-4a 

Sulfuric  Acid  (H2SO4)  Neutron-Diffraction  10°K  Geometry 

NDb 

Spc 

SEC 

MF 

Ipc 

S-Ol 

1.48 

1.554 

1.551 

1.549 

1.552 

S=02 

1.49 

1.425 

1.428 

1.431 

1.428 

Ol-H 

0.91 

0.972 

0.976 

0.980 

0.979 

01...02d 

2.52 

2.55 

2.55 

2.55 

2.55 

H...02d 

1.65 

1.59 

1.58 

1.58 

1.58 

ZH-Ol-S 

113.6 

113.5 

114.4 

114.3 

115.2 

zoi-s-or 

112.8 

104.9 

105.7 

106.0 

105.4 

Z02-S-02' 

110.1 

120.0 

119.1 

118.6 

119.5 

ZOl-S-02 

108.0 

107.2 

107.2 

107.0 

107.3 

Z01-S-02' 

109.0 

108.3 

108.4 

108.8 

108.2 

Z01-H...02d 

175 

171.5 

173.8 

171.3 

172.7 

ZS-01...02d 

111.9 

112.1 

112.0 

112.1 

112.3 

Z01...02d-Sd 

155.3 

159.3 

159.1 

159.1 

158.9 

x(H-Ol-S-Or) 

84.1 

93.1 

90.5 

92.6 

91.3 

T(H-01-S-02) 

-36.5 

-21.9 

-25.1 

-23.4 

-23.9 

Energy: 

Hartree 

-698.136173 

-698.165534 

kcal/mole 


-698.160937  -698.160931 

-438,086.963  -438,105.387 

-438,102.503  -438,102.499 
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a.  Bond  lengths  in  angstroms.  Bond  angles  in  degrees. 

b.  Neutron  diffraction  10°K  structure  from  reference  37.  The  estimated 
uncertainties  are  reported  in  the  paper  as  ±0.05-0.10  A. 

c.  This  work.  Calculated  using  the  6-3 IG**  basis  set. 

d.  Atom  is  located  at  fractional  coordinates  -l/2+x,l/2-i-y,z  on  an  adjacent 
molecule. 


47 


model,  the  remaining  differences  between  the  calculated  and  experimental  crystal 
structures  need  closer  scrutiny. 

The  differences  between  the  calculated  and  neutron-diffraction  crystal 
structures  are  easily  explained  by  examining  the  experimental  data.  Moodenbaugh 
and  coworkers  report  estimated  standard  deviations  (esd's)  of  ±  0.05-0.10  A. 
Although  their  calculated  esd's  are  about  0.01  A,  the  larger  esd's  arise  from  an 
apparent  discrepancy  in  the  01 -H  bond  length  which  is  known  to  be  temperature 
independent  to  about  ±  0.002  A. 37  Based  on  their  ±  0.05-0. 10  A  uncertainty  in 
bond  lengths,  all  of  the  calculated  values  fall  within  that  range.  Any  comparison 
with  the  neutron-diffraction  study  by  Moodenbaugh  and  coworkers  is  essentially 
useless  because  of  the  extremely  large  uncertainties  in  the  experimental  data. 

There  are  numerous  reasons  for  the  large  errors  reported  by  Moodenbaugh 
and  coworkers.  First,  the  structural  model  that  is  used  in  any  refinement  of  the  raw 
crystal  data  is  crucial  to  the  results.  In  their  study,  they  based  the  heavy  atom 
positions  on  the  1965  X-ray  study  with  its  larger  esd's  instead  of  the  1978  study 
with  its  smaller  esd's.  They  were  unaware  of  the  1978  study.  They  also  chose  to 
use  isotropic  temperature  factors  for  simpUcity,  even  though  both  of  the  X-ray 
diffraction  studies  indicate  that  anisotropic  factors  would  be  more  appropriate. 
Their  R  factors,  on  the  order  of  0.20,  are  approximately  double  those  from  the  X- 
ray  diffraction  studies.  Second,  the  quaUty  of  their  data  is  rather  poor  with  an 
unfavorable  typical  peak  to  background  ratio  of  1:2.  This  is  probably  due  to  the 
fact  that  they  used  a  powder  sample  since  they  could  not  obtain  a  suitable  single 
crystal.  The  high  background  caused  by  incoherent  scattering  required  them  to  use 
the  Rietveld  refinement  method  which  is  not  normally  used. 
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Finally,  the  data  from  the  neutron-diffraction  experiment  is  questionable  at 
best.  Table  3-5  shows  the  bond  lengths  and  angles  reported  by  them  at  240°K. 
However,  when  I  used  their  coordinates  to  calculate  structural  parameters,  the 
actual  S-Ol  bond  length  differs  greatly  from  that  they  reported.  At  first  glance, 
only  the  reported  S-01  bond  length  would  appear  to  be  in  error.  If  the  value  of 
1.67  A  that  I  calculated  is  correct,  the  S-Ol  bond  length  would  undergo  a 
shortening  of  0.19  A  between  240°K  and  10°K  which  is  highly  unlikely  since 
H2SO4  was  found  to  be  ordered  at  both  temperatures.  In  addition,  the  values  I 
calculate  from  their  data  for  the  ZOl-S-01'  (=87.9°)  and  Z02-S-02'  (=133.1°) 
differ  gready  from  those  in  the  X-ray  diffraction  study  (103.7°  and  1 18.3°)  and  in 
the  neutron-diffraction  study  at  10°K  (1 12.8°  and  1 10.1°).  The  discrepancies  in  the 
S-Ol  bond  length  and  the  angles  about  sulfur  suggest  to  me  that  there  is  an  error  in 
the  fractional  coordinates  of  the  sulfur  atom.  Since  sulfur  occupies  the  special 
position  0,y,  1/4,  the  error  appears  to  be  in  the  y  fractional  coordinate.  I  have  made 
several  attempts  to  determine  the  correct  sulfur  coordinates  by  varying  the  y 
coordinate.  In  my  opinion,  the  most  probable  structure  for  H2SO4,  based  on  the 
data  presented  in  the  paper  at  240°K,  occurs  when  y =0.068  as  shown  in  Table  3-5. 
The  reported  value  is  y=0.028(6).  My  value  of  0.068  appears  to  be  a  logical 
transposition  of  the  reported  value.  At  y=0.068,  there  is  much  better  agreement 
between  the  10°K  and  240°K  structures.  I  have  talked  to  Moodenbaugh  about  this 
problem.  Moodenbaugh  was  not  previously  aware  of  the  problem  and  told  me  he 
put  much  more  rehance  on  the  heavy  atom  positions  from  the  X-ray  studies. 
Further,  he  indicated  absolutely  no  interest  in  trying  to  resolve  the  problem.  The 
numerous  problems  with  the  neutron-diffraction  study  by  Moodenbaugh  and 
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Table  3-5^ 

H2SO4  Neutron-Diffraction  Structural  Parameters  at  240°K 


Reported'' 

Calculated 

Probable 

S-Ol 

1.54 

1.67 

1.54 

S=02 

1.40 

1.40 

1.49 

Ol-H 

1.00 

1.00 

1.00 

01... 02c 

2.47 

2.47 

2.47 

H...02C 

1.49 

1.50 

1.50 

Z01-H...02C 

163 

164 

164 

ZOl-S-Or  87.9  97.9 

Z02-S-02'  133.1  119.6 

a.  Bond  lengths  in  angstroms.  Bond  angles  in  degrees. 

b.  Neutron  diffraction  10°K  structure  from  reference  37.  The  estimated 
uncertainties  are  reported  in  the  paper  as  ±0.05-0.10  A. 

c.  Atom  is  located  at  fractional  coordinates  -l/2+x,l/2-i-y,z  on  an  adjacent 
molecule. 


50 


coworkers  certainly  warrant  another  neutron-diffraction  study  of  H2SO4. 

A  comparison  of  the  calculated  and  X-ray  crystal  structures  in  Table  3-3 
shows  relatively  good  agreement.  By  far  the  largest  discrepancy  is  in  the  S-Ol 
single  bond  length.  The  X-ray  experiment  determined  the  S-Ol  distance  to  be 
about  0.03  A  shorter  than  the  calculation.  As  discussed  earlier,  a  major 
disadvantage  of  X-ray  diffraction  is  that  it  determines  distances  between  centers  of 
electron  density  and  not  nuclei.  Since  the  shorter  bond  length  could  be  a 
consequence  of  the  X-ray  experiment,  I  decided  to  investigate  it  further.  In  order  to 
investigate  the  S-Ol  bond  length  discrepancy,  I  chose  three  distinct  approaches. 
The  first  was  to  calculate  the  crystal  structure  with  an  expanded  basis  set.  The 
second  involved  calculations  on  the  H2SO4  dimer.  In  the  third  method,  I  used  the 
program  MOCALC^^  in  order  to  look  at  the  total  electron  density  along  the  S-Ol 
bond. 

Since  the  point  charge  model  may  not  adequately  describe  the  hydrogen 
bonding,  the  calculated  S-Ol  bond  length  may  be  too  long.  However,  various 
studies  on  dimers  have  shown  that  not  only  similar  double  zeta  basis  sets'^9,5 1,54,55 
but  also  even  the  minimal  STO-3G  basis  set^^  adequately  describe  hydrogen 
bonding.  Saeb0  and  coworkers  concluded  from  their  study  of  cyanoformamide 
that  hydrogen  bonding  was  dominated  by  electrostatic  effects.^  In  addition,  Tang 
and  Fu  concluded  from  ab  initio  studies  of  hydrogen  bonding  between  methyl 
cyanide  (CH3CN)  or  methyl  isocyanide  (CH3NC)  and  methanol  (CH3OH)  that  it 
was  a  result  of  ionic  binding  of  two  closed- shell  systems.'^^  Nevertheless,  in  order 
to  see  if  an  expansion  in  the  basis  set  might  better  describe  the  hydrogen  bonding,  a 
calculation  was  made  with  the  6-31  -i-G**  basis  set.^O  The  "+"  implies  the  addition 
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of  a  diffuse  sp  shell  to  the  hydrogen  atoms  (exponent=0.0360).  The  addition  of 
diffuse  functions  is  required  in  the  calculation  of  anions,  transition  structures,  and 
other  species  with  significant  electron  density  far  removed  from  the  nuclear  center. 
The  optimized  geometry  with  the  6-31  -i-G**  basis  set  was  practically  identical  to 
that  of  the  6-3 IG**  basis  set  with  no  difference  in  the  S-Ol  bond  length. 

Another  approach  I  pursued  was  to  optimize  the  structure  of  the  H2SO4 
dimer.  The  two  molecules  chosen  were  those  that  were  hydrogen  bonded  in  the 
crystal.  The  dimer  calculation  showed  a  definite  shortening  of  the  S-Ol  bond  at  the 
STO-3G*  level.  The  hydrogen  bonded  S-Ol  and  S-02  bond  lengths  were  1.58  A 
and  1.46  A,  respectively,  compared  to  1.62  A  and  1.45  A  in  the  monomer  at  the 
ST0-3G*  level.  The  shortening  shown  in  the  dimer  must  be  viewed  with  caution. 
In  the  optimized  dimer,  the  two  molecules  rotated,  as  a  result  of  the  program,  such 
that  two  hydrogen  bonds  were  formed  between  them  instead  of  only  one  as  found 
in  the  crystal.  This  is  apparently  caused  by  the  lack  of  neighbor  molecules  which 
prevent  the  rotation  in  the  crystal.  My  attempts  to  conduct  dimer  calculations  at  the 
6-3 IG**  level  with  the  inclusion  of  the  point  charges  have  not  been  successful.  I 
initially  planned  to  optimize  the  dimer  structure  at  the  ST0-3G*  level  with  the 
inclusion  of  the  point  charges.  However,  because  of  the  success  of  the  total 
electron  density  plot  in  explaining  the  shortened  bond  in  the  X-ray  study,  I  have  not 
completed  any  additional  dimer  calculations. 

The  program  MOCALC^^  allows  plotting  electron  density  of  either  a  single 
molecular  orbital  or  all  of  the  molecular  orbitals.  MOCALC  was  specifically 
written  to  use  the  data  generated  in  the  TEXAS  program.  In  order  to  investigate  the 
shortened  S-Ol  bond,  I  have  completed  a  total  electron  density  plot  of  H2SO4  in 
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the  Ol-S-Or  plane  which  is  shown  in  Figure  3-3.  As  Figure  3-3  clearly  shows, 
the  electron  density  of  the  Ol  atom  is  definitely  shifted  toward  the  S  atom.  After 
talking  with  several  crystallographers,  I  have  determined  that  there  is  no  commonly 
used  quantitative  method  to  account  for  this  apparent  bond  shortening.  The 
crystallographers  simply  accept  the  fact  that  it  exists.  In  order  to  get  a  quantitative 
feel  for  the  magnitude  of  this  shortening,  I  have  simply  assumed  the  location  of  the 
center  of  the  electron  density  on  the  Ol  atom  in  Figure  3-3.  By  making  simple 
measurements  with  a  ruler,  I  have  determined,  that  in  the  case  of  H2SO4,  the 
shortening  of  the  S-Ol  bond  length  in  the  X-ray  study  is  on  the  order  of  0.03-0.05 
A.  Thus,  as  shown  in  Figure  3-3, 1  conclude  that  the  short  S-Ol  bond  length 
reported  in  the  X-ray  studies  is  a  consequence  of  the  X-ray  experiment. 

The  various  methods  used  in  this  study  to  determine  the  magnitude  and 
location  of  the  point  charges  yielded  some  interesting  results.  As  shown  in  Table 
3-3  and  3-4,  the  SF  and  SE  methods  yielded  almost  identical  geometries  despite  the 
differences  in  the  magnitude  of  the  point  charges.  In  both  cases,  the  magnitude  of 
the  point  charges  obtained  from  the  experimental  crystal  structures  in  the  SE 
method  differed  by  as  much  as  15%  from  those  calculated  for  the  optimized  free 
molecule.  The  differences  between  the  SF  and  SE  structures  were  no  more  than 
0.004  A  in  bond  lengths  and  0.9°  in  bond  angles.  In  most  cases,  within  the 
uncertainties  of  the  calculations,  the  two  methods  gave  structures  which  agree.  The 
convincing  agreement  between  these  two  methods  shows  that  the  structures 
calculated  using  the  point  charge  method  are  relatively  insensitive  to  the  magnitude 
of  the  charges.  This  indicates  that  the  arbitrary  manner  in  which  the  MuUiken 
population  analysis  assigns  net  charge  has  httle  effect  on  the  calculated  results. 
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Figure  3-3 

Total  Electron  Density  Plot  of  H2SO4  in  the  01-S-Ol*  Plane 


CONTOUR    INTERVAL    -0.060000 
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My  study  shows  that  the  iterative  technique  used  in  the  MF  and  IF  methods 
offers  several  advantages  over  the  single  optimization  cycle  methods.  Both  of  these 
methods  show  substantial  improvements  in  the  total  energy  over  the  SF  method. 
The  energy  improvements  are  on  the  order  of  15-19  kcal/mole.  Throughout  the  MF 
and  IF  optimization  cycles,  the  calculated  structural  parameters  continued  to  better 
approximate  those  in  the  crystal.  Specifically,  the  S-Ol  bond  length  continued  to 
shorten,  the  ZOl-S-Ol'  increased,  and  the  Z02-S-02'  decreased.  The 
improvements  in  both  the  energy  and  structural  parameters  confirm  the  use  of  the 
iterative  approach.  One  particularly  interesting  result  from  the  MF  and  IF  methods 
is  the  change  in  gross  atomic  populations.  Three  ab  initio  studies  by  Ottersen  and 
Jensen  found  that  there  is  a  transfer  of  electron  density  from  the  hydrogen  atom  to 
the  acceptor  atom.  The  formamide  (HCONH2)  dimer  showed  a  transfer  of  0.05 
electron  across  the  hydrogen  bond.54  Ottersen's  study  of  the  keto  (HCONH2)  and 
enol  (HOCHNH)  tautomers  of  formamide  exhibited  gross  atomic  population 
increases  of  0.057  and  0.089  on  the  acceptor  atoms  and  decreases  of  0.046  and 
0.044,  respectively,  in  the  bridging  hydrogens.^^  Their  study  of  the  four 
monohydrated  formamide  complexes  demonstrated  gross  atomic  population 
increases  of  0.009-0.034  for  the  acceptor  atoms  and  decreases  of  0.015-0.039  for 
the  bridging  hydrogens.^l  The  gross  atomic  populations  for  my  calculations  using 
the  iterative  method  are  presented  in  Table  3-6.  All  four  of  my  calculations  clearly 
follow  the  trend  of  Ottersen  and  Jensen  by  exhibiting  a  transfer  of  electron  density 
from  the  bridging  hydrogens  to  the  acceptor  atoms  (02/02').  These  results 
confirm  the  abiUty  of  the  point  charge  model  to  accurately  describe  the  crystal 
environment  for  hydrogen  bonded  systems.  In  summary,  the  iterative  technique 
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Table  3-6a 

H2SO4  Gross  Atomic  Populations 


Free 

Optimized 

Net 

Molecule 

Geometry^ 

Change 

X-ray: 

MF-               S 

14.18843 

14.11860 

-0.070 

01/or 

8.65419 

8.67702 

-hO.023 

H/H' 

0.60258 

0.51039 

-0.092 

02/02' 

8.64902 

8.75329 

+0.104 

IF-                 S 

14.18843 

14.10879 

-0.080 

01/or 

8.65419 

8.67743 

-hO.023 

H/H' 

0.60258 

0.50813 

-0.094 

02/02' 

8.64902 

8.76004 

+0.111 

Neutron-Diffraction: 

MF-               S 

14.18843 

14.11909 

-0.069 

01/01' 

8.65419 

8.66920 

+0.015 

H/H' 

0.60258 

0.50504 

-0.098 

02/02' 

8.64902 

8.76621 

+0.117 

IF- 


S 

01/01' 

H/H' 

02/02' 


14.18843 
8.65419 
0.60258 
8.64902 


14.10930 
8.67745 
0.50761 
8.76029 


-0.079 
+0.023 
-0.095 
+0.111 


a.  All  values  were  calculated  from  the  Mulliken  population  analysis  with  the  6- 

31G**  basis  set. 

b.  These  values  were  obtained  using  the  indicated  methods  after  the  calculations 

converged  in  both  energy  and  structure.  In  all  cases  this  required  5 
optimization  cycles. 
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tends  to  minimize  any  detrimental  effects  of  the  Mulliken  population  analysis.  In 
addition,  by  allowing  the  point  charge  model  to  dynamically  change  in  relation  to  its 
environment,  any  effects  or  forces  which  are  not  entirely  electrostatic  by  nature 
appear  to  be  accounted  for.  The  only  disadvantage  of  the  iterative  technique  is  the 
amount  of  computation  time  required.  With  each  run  averaging  2,500  seconds  and 
up  to  30  runs  required,  the  iterative  technique  is  fairly  expensive.  However,  the 
accurate  results  from  these  methods  far  outweigh  their  cost.  The  MF  and  IF 
methods  are  clearly  superior  to  the  SF  method  used  by  Saeb0  and  coworkers. 

My  study  has  shown  that  the  IF  method  is  particularly  useful.  Not  only 
does  it  allow  calculations  on  molecules  for  which  incomplete  experimental 
structures  exist,  but  it  also  gave  better  results  in  my  calculations  based  on  the  X-ray 
data.  The  fact  that  the  MF  and  IF  X-ray  structures  are  almost  identical  substantiates 
the  point  charge  model.  The  point  charge  model  reproduces  the  deformation  forces 
which  are  present  in  the  solid  upon  crystallization  as  shown  by  the  IF  method.  In 
the  IF  method,  the  point  charges  are  initially  positioned  in  the  free  molecule 
geometry  but  they  are  subsequendy  adjusted  to  the  newly  calculated  "free"  molecule 
geometry.  Since  the  point  charges  eventually  assume  the  approximate  structure  of 
the  crystal,  the  point  charge  model  must  essentially  reproduce  the  forces  present  in 
the  crystal.  Based  on  my  results,  both  the  MF  and  IF  methods  are  superior  to  the 
SF  method  and  they  allow  much  more  flexibility  in  the  model.  Since  both  the  MF 
and  IF  methods  have  been  shown  to  give  almost  identical  results,  either  method  can 
be  used  but  the  IF  method  should  certainly  be  used  in  cases  where  the  experimental 
data  is  incomplete. 

A  major  concern  with  the  point  charge  model  is  the  dependence  of  the 
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calculated  structure  on  the  magnitude  of  the  point  charges.  In  order  to  investigate 
the  dependance,  I  made  several  runs  in  which  I  varied  only  the  magnitude  of  the 
point  charges.  The  point  charges  were  positioned  according  to  the  arrangement 
found  in  the  10°K  neutron-diffraction  study  and  the  free  molecule  was  optimized 
inside  these  point  charges  as  in  the  SF  method.  The  results  of  those  calculations 
are  shown  in  Table  3-7.  In  run  #1  which  uses  the  SE  method,  the  magnitude  of  the 
point  charges  was  calculated  from  the  experimental  geometry.  The  net  charges  in 
run  #2  were  calculated  from  the  optimized  free  molecule.  Run  #2  is  equivalent  to 
the  SF  results  reported  in  Table  3-4.  In  run  #3, 1  arbitrarily  normalized  the  net 
charges  of  the  free  molecule  to  the  02/02'  net  charge.  Despite  the  difference  in  the 
net  charges  which  vary  by  up  to  80%,  the  calculated  geometries  are  extremely 
similar.  The  maximum  difference  in  bond  lengths  and  bond  angles  is  0.02  A  and 
2.1°,  respectively,  occurring  only  in  run  #3  with  its  extreme  values  of  charge.  I 
conclude  from  these  results  that  the  structures  calculated  using  the  point  charge 
model  are  fairly  insensitive  to  the  magnitude  of  the  point  charges.  Therefore,  the 
arbitrary  assignment  of  net  charge  by  the  Mulliken  population  analysis  has  little 
effect  of  the  calculated  structure. 
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Table  3-7^ 
Calculated  H2SO4  Solid  Phase  Geometries 

#lb  #2b  #3b 


Net  Charge: 

S 

1.75069 

1.81157 

2.79124 

oi/or 

-0.56293 

-0.65419 

-1.00796 

H/H- 

0.39592 

0.39742 

0.61234 

02/02' 

-0.70833 

-0.64902 

-1.00000 

Structure: 

S-Ol 

1.551 

1.554 

1.549 

S=02 

1.428 

1.425 

1.436 

Ol-H 

0.976 

0.972 

0.992 

ZH-Ol-S 

114.4 

113.5 

115.6 

zoi-s-or 

105.7 

104.9 

106.6 

Z02-S-02' 

119.1 

120.0 

118.1 

ZOl-S-02 

107.2 

107.2 

106.4 

ZOl-S-02' 

108.4 

108.3 

109.4 

X(H-Ol-S-Ol') 

90.5 

93.1 

91.1 

T(H-01-S-02) 

-25.1 

-21.9 

-25.6 

Energy: 

Hartree 

-698.160937 

-698.136173 

-698.194785 

kcal/mole 

-438,102.503 

-438,086.963 

-438,123.743 
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a.  Bond  lengths  in  angstroms.  Bond  angles  in  degrees. 

b.  Geometry  optimized  at  the  6-3 IG**  level.  The  point  charge  locations  are  based 
on  the  10°K  neutron-diffraction  data.  The  magnitude  of  the  point  charges  were 
obtained  as  explained  in  the  text. 


Chapter  Four 

CONCLUSION 

In  this  paper,  I  have  presented  a  thorough  ab  initio  study  of  sulfuric  acid 
and  the  point  charge  model.  The  point  charge  model  has  been  shown  to  be 
extremely  advantageous  in  the  calculation  of  crystal  structures  with  few  limitations. 
In  addition,  I  have  presented  a  comprehensive  review  of  various  crystal  field 
simulation  methods  which  have  been  used.  In  this  chapter,  the  important  results 
from  my  study  are  summarized  and  my  suggestions  for  future  work  are  offered. 

My  calculations  of  the  gas  phase  structure  of  H2SO4  show  that  the  6-3 IG** 
basis  set  accurately  reproduces  the  microwave  structure.  With  few  exceptions,  the 
calculated  6-31G**  structure  is  within  the  experimental  uncertainties  of  the 
microwave  study.  The  4-31G**  basis  set  gave  a  similar  structure  to  the  6-31G** 
results  with  a  substantial  savings  in  computation  time.  My  study  has  confirmed 
that  the  use  of  displaced  p  functions  to  simulate  d  functions  has  essentially  no  effect 
on  the  calculated  structure.  The  only  difference  being  in  the  total  energy  which  is 
less  than  1  kcal/mole  lower  with  the  true  d  functions.  The  ab  initio  programs 
TEXAS(D)  and  Gaussian-82  give  identical  results  with  the  6-3 IG**  basis  set  for 
H2SO4  as  expected.  My  calculations  confirm  that  polarization  functions  are 
absolutely  essential  to  reproduce  torsional  angles  on  oxygen. 

My  calculations  of  the  crystal  structure  of  H2SO4  demonstrate  that  the  point 
charge  model  accurately  reproduces  the  structural  trends  which  occur  in 
transforming  from  the  gas  to  soHd  phase.  My  study  of  the  point  charge  model 
proves  that  it  precisely  simulates  the  defonnation  forces  which  are  present  in  the 
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solid  upon  crystallization.  During  the  course  of  my  study,  I  discovered  that  the 
only  neutron-diffraction  study  of  H2SO4  is  plagued  by  numerous  errors,  rendering 
it  totally  useless. 

The  calculated  crystal  structures  based  on  the  X-ray  diffraction  study  are  in 
good  agreement  with  the  experiment.  The  major  discrepancy  is  in  the  S-O  single 
bond  length  which  is  calculated  about  0.03  A  longer  than  in  the  X-ray  study.  After 
a  thorough  investigation,  I  conclude  that  the  discrepancy  is  due  to  a  shifting  of 
electron  density  from  the  oxygen  atom  towards  the  sulfur  atom,  resulting  in  the 
short  experimental  bond  length. 

The  iterative  technique  gives  superior  results  to  any  of  the  single 
optimization  cycle  methods.  Even  with  the  point  charges  arranged  in  the  gas 
geometry,  the  point  charge  model  gives  almost  identical  structural  results  to  those 
calculations  where  the  charges  are  positioned  in  the  crystal  structure.  The  point 
charge  model  more  than  adequately  describes  the  crystal  forces,  especially  for 
systems  involving  hydrogen  bonding. 

The  iterative  technique  minimizes  any  possible  drawbacks  of  the  point 
charge  model.  Any  consequence  of  the  arbitrary  assignment  of  net  charge  by  the 
Mulliken  population  analysis  is  minimized  by  allowing  the  charges  to  vary.  By 
allowing  the  point  charge  model  to  dynamically  change  in  relation  to  its 
environment,  any  forces  which  are  not  electrostatic  by  nature  appear  to  be 
accounted  for. 

Further  evidence  of  the  point  charge  model's  ability  to  represent  the  crystal 
environment  is  exhibited  in  the  change  in  gross  atomic  population  throughout  the 
iterative  cycles.  There  is  a  definite  shift  of  electron  density  from  the  bridging 
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hydrogens  to  the  acceptor  atoms  (02/02'),  identical  to  those  found  in  other  ab  initio 
studies  of  dimers. 

My  study  proves  that  the  calculated  crystal  structures  are  insensitive  to  the 
magnitude  of  the  point  charges.  I  conclude  that  the  arbitrary  assignment  of  charge 
by  the  Mulliken  population  analysis  has  a  negligible  effect  on  the  point  charge 
model. 

SUGGESTIONS  FOR  FUTURE  WORK 

The  remainder  of  this  chapter  is  a  short  presentation  of  areas  that  I  feel  need 
further  work. 

Although  the  results  of  my  study  confirm  the  validity  of  the  point  charge 
model,  the  only  types  of  crystals  that  have  been  tested  to  date  involve  hydrogen 
bonding.  I  strongly  feel  that  further  studies  should  be  done  on  crystals  bound 
together  by  both  ionic  and  covalent  bonds. 

A  quantitative  method  needs  to  be  developed  to  treat  the  shift  of  electron 
density  that  frequently  occurs  in  X-ray  diffraction  studies.  Currendy,  no 
quantitative  method  exists  and  bond  lengths  determined  by  X-ray  diffraction  are 
often  reported  shorter  than  their  actual  value.  The  development  of  such  a  method  is 
desperately  needed  to  facilitate  comparisons  between  calculated  and  experimental 
crystal  structures. 

The  problems  that  I  have  presented  with  the  neutron-diffraction  study  of 
H2SO4  by  Moodenbaugh  and  coworkers  point  to  repeating  that  study.  The  new 
study  should  incorporate  the  Yu  and  Mak  data  and  be  done  on  a  single  crystal  if 
possible.  As  a  minimum,  Moodenbaugh's  data  should  be  refined  again  using  the 
heavy  atom  positions  from  the  1978  X-ray  study  and  anisotropic  temperature 
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factors.  I  expect  that  the  results  of  a  new  neutron-diffraction  study  will  confirm 
that  the  S-O  single  bond  length  is  closer  to  my  calculated  value  than  the  X-ray 
value. 

After  discussions  with  several  crystallographers,  they  all  feel  that  a  present 
day  X-ray  diffraction  study  of  H2SO4  could  locate  the  hydrogens.  With  the  aid  of 
a  crystallographer,  I  have  unsuccessfully  tried  to  refine  the  Yu  and  Mak  data  to 
locate  the  hydrogens.  Since  no  further  knowledge  can  be  gained  from  past  X-ray 
diffraction  studies,  I  feel  a  new  X-ray  study  should  be  completed.  The  positions  of 
the  hydrogens,  if  located,  can  then  be  compared  to  the  neutron-diffraction  study  to 
establish  its  accuracy. 
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